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ABSTRACT 


This  thesis  is  a  study  of  missile  and  target  parameters 
used  in  second  and  third  order  modeling  of  the  tracking 
subsystem  used  in  radar  guided  missiles.  Guidance  methods 

are  analyzed  to  determine  which  method  is  optimum  in  a 

<r  c  '■ 

search  for  an /(nidealw/ missile .  Target  parameters  which  have 
an  effect  on  the  missile  tracking  system  are  analyzed  and  a 
target  acceleration  probability  model  is  discussed.  A  two 
dimensional  third  order  tracking  model  is  simulated 
utilizing  a  Kalman  filter  for  target  parameter  estimation 
and  prediction.  Linear  second  and  third  order  tracking 
models  are  simulated  and  compared  with  the  third  order 
Kalman  filter  tracker.  This  thesis  concludes  that  a 
proportional  navigation  guidance  method,  with  a  non  linear 
third  order  tracking  Kalman  filter,  is  the  better  model. 
Benefits  of  using  a  non  linear  third  order  Kalman  filter  may 
not  overide  the  cost  and  complexity  of  implementation  of  the 
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I.  INTRODUCTION 

Aviation  plays  an  extensive  role  in  current  combat 
scenarios.  An  aircraft,  because  of  its  capability  to  carry 
missiles,  is  a  very  formidable  weapon  platform.  Missiles 
provide  offensive  killing  power  which  change  tactics  in 
battle  scenarios.  In  order  to  have  the  edge  in  the  air  to 
air  arena,  an  aircraft  must  possess  the  best  defensive  and 
offensive  capabilities.  One  of  the  main  weapons  in  the 
aircraft  arsenal  is  the  radar  guided  missile. 

Radar  guided  missiles  are  all  weather  capable,  can  be 
launched  outside  of  visual  range  and  are  less  susceptible  to 
countermeasures,  compared  with  other  missile  types. 
Improvement  of  the  missile  is  a  constant  necessity  to 
maintain  air  superiority. 

Improvements  in  aircraft  maneuverability  dictate  the 
need  for  missiles  to  increase  performance  and  capabilities. 
A  rule  of  thumb  for  design  is  that  the  missile  must  have  a 
4:1  acceleration  advantage  over  the  target.  With  modern 
aircraft  able  to  sustain  lateral  accelerations  of  ten  times 
the  force  of  gravity  (G)  the  missile  must  be  capable  of  40G. 

Guidance  methods  are  chosen  which  optimize  the  missile 
capabilities  to  destroy  the  target.  The  guidance  method 
selected  has  a  large  impact  on  the  design  of  other  missile 
subsystems.  For  any  missile  to  guide  to  the  target,  the 
sensor  subsystem  must  track  the  target.  To  optimize  the 
guidance  and  tracking  of  radar  guided  missiles  a  predicting 
filter  can  be  used. 

One  of  the  simplest  missile  guidance  techniques  is  a 
second  order  system  which  compensates  for  bearing  error  and 
bearing  rate  error.  This  thesis  will  look  at  third  order 
models  to  help  optimize  the  missile  sensor  subsystem  to 
provide  better  guidance  command  inputs.  A  major  impetus  for 
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finding  an  optimum  predicting  guidance  method  is  to  improve 
missile  performance  in  the  final  guidance  stages. 

There  is  a  region  at  the  end  of  the  missile  flight  path, 
where  the  time  to  intercept  is  so  short  that  inputs  to  the 
missile  control  surfaces  will  not  be  effective.  If  the 
missile  has  an  exact  solution  of  target  parameters,  it  can 
predict  the  future  target  position,  through  the  time  where 
no  inputs  will  be  effective.  The  missile  is  flying  to  the 
projected  position  of  the  target,  compensating  for  the  time 
delay  of  control  effectiveness. 

Section  II  will  look  at  some  guidance  methods  and  target 
parameters  to  aid  in  finding  the  optimum  missile  guidance. 
Section  III  looks  at  target  models  and  how  to  implement  them 
in  the  computer  simulation.  Section  IV  derives  the  second 
and  third  order  models.  A  Kalman  Predicting  Filter  is 
discussed  in  Section  V.  Computer  simulation  and 
implementation  of  two  dimensional  third  order  models  with 
the  Kalman  Filter  is  given  in  Section  VI.  Section  VII 
contains  conclusions  and  recommendations.  Program  listings 
of  the  computer  simulation  are  in  Appendix  A  and  Appendix  B. 
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II.  THE  IDEAL  MISSILE 

Definition  of  an  ideal  missile  is  very  difficult.  Cost 
or  performance  functions  can  be  generated  to  account  for 
miss  distance,  fuel  expended,  control  inputs,  flight  path 
and  numerous  other  parameters.  To  the  operator  the  ideal 
missile  is  one  that  destroys  the  target.  The  designer  must 
try  to  account  for  a  variety  of  scenarios  and  targets  to 
design  the  optimum  missile.  Tradeoffs  of  performance  and 
costs  will  dictate  what  the  final  sub-optimum  missile  will 
be . 

In  an  attempt  to  find  an  ideal  missile  insight  may  be 
gained  by  looking  at  the  various  flight  paths  and  parameters 
for  different  guidance  methods.  Three  methods  to  look  at 
are  pure  pursuit,  lead  pursuit  and  proportional  navigation. 
Pure  pursuit  would  entail  a  missile  always  flying  directly 
at  the  target.  Lead  pursuit  would  be  the  case  of  a  missile 
always  flying  to  a  point  slightly  ahead  of  the  target.  The 
magnitude  of  the  lead  may  vary  from  one  time  constant,  to 
the  total  missile  flight  time,  ahead  of  the  target.  The 
latter  would  produce  a  lead  collision  intercept. 
Proportional  navigation  utilizes  line  of  sight  rate  to  guide 
the  turn  of  the  missile,  to  zero  any  further  line  of  sight 
rate.  The  commanded  acceleration  to  the  missile  is  given  by 
the  equation: 

AM  =  k  Me  (1) 

where  AM  =  missile  acceleration 

k  *  constant  of  proportionality 
£  =  antenna  angle  rate 

Vc  *  closing  velocity 


The  constant  of  proportionality  is  determined  by  the 
designer.  Observing  the  effects  of  the  proportionality 


constant  on  the  line  of  sight  rate  shown  in  Figure  II-l,  any 
constant  above  3  is  an  appropriate  value.  For  k  less  than  3 
the  line  of  sight  rate  has  its  largest  slope  at  the  end  of 
the  intercept.  For  k  greater  than  3  the  line  of  sight  rate 
will  be  very  small  prior  to  the  impact  point. 


Figure  II-l.  Frop  Nav  Proportionality  Constant  [Ref.  1] 

Three  guidance  methods  are  analyzed  to  compare  flight 
paths,  heading  changes,  line  of  sight  angle  and  line  of 
sight  rate.  A  fourth  missile  is  simulated  and  plotted, 
identified  by  "direct". 

The  direct  path  missile  is  programmed  using  "a 
posteriori"  knowledge  of  the  target  flight  path.  The  direct 
path  missile  goes  to  the  point  of  impact,  in  a  straight 
line,  from  the  point  of  launch. 

If  the  "ideal  missile"  is  defined  such  that  there  is  no 
missile  maneuvering,  burns  minimal  fuel,  has  the  greatest 
launch  distance,  minimum  intercept  time  and  accounts  for  all 
target  maneuver,  it  will  be  a  direct  path  missile. 


Two  programs  were  written  using  Dynamic  Simulation 

Language,  DSL,  to  analyze  the  guidance  methods  and  produce 
plots.  The  programs  are  included  in  Appendix  A. 

Three  scenarios  are  used  to  compare  the  guidance 
methods : 

head-on  aspect 
tail  aspect 
beam  aspect 

Three  target  accelerations  are  used  to  determine  the 
effect  of  the  target  maneuver  on  the  parameters.  The 

selected  target  accelerations  are: 

0  g '  s 
3  g '  s 
6  g '  s 

A  second  order  proportional  navigation  missile  with  a 
proportional  navigation  guidance  constant  of  four,  where 
the  line  of  sight  angle  and  angle  rate  is  estimated  by  the 
antenna  parameters  of  angular  position  and  angular  rate. 
The  governing  equations  for  tracking  in  bearing  are: 

3  =  p  dt  +  p„  (2) 

P  =  (T  dt  +  p0  (3) 

P’=  -20*B  +  100* (LOS-p )  (4) 

where  P  *  antenna  position  angle 
P  =  antenna  angle  rate 

P  =  antenna  angular  acceleration 

LOS  «  actual  target  bearing 

Pure  pursuit  and  lead  pursuit  guidance  missiles  are 
initialized  heading  directly  at  the  target.  Pure  pursuit 
guidance  maintains  heading  directly  at  the  target  by 
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calculating  the  heading  at  each  step  of  the  discrete 
simulation.  The  equation  for  pure  pursuit  heading  is: 

PPHDG  *  atan ( (yt-ym) / <xt-xm) )  (5) 

where  yt  *  current  target  Y  coordinate 
xt  *  current  target  X  coordinate 
ym  =  current  missile  Y  coordinate 
xm  *=  current  missile  X  coordinate 

Lead  pursuit  maintains  a  heading  in  front  of  the  target, 
with  a  variable  lead,  calculated  using  half  the  time  to  go, 
given  by: 

LPHDG  =  atan ( (yt+tvely* . 5ttg) / (xt+tvelx* . 5*ttg) )  (6) 

where  yt  =  current  target  Y  coordinate 

xt  =  current  target  X  coordinate 

tvely  =  target  velocity  in  Y  direction 
tvelx  =  target  velocity  in  X  direction 
ttg  *  time  to  go  in  the  intercept 

The  results  of  the  comparisons  are  given  in  graph  form 

in  Figure  II-2  through  Figure  11-37. 

A.  HEAD-ON  ASPECT 

Figures  II-2  through  11-13  are  the  results  of  missile 
guidance  comparisons  for  head-on  aspect  initial  condition 
with  OG,  3G  and  6G  constant  target  acceleration.  The 
missile  begins  at  the  origin  of  the  graph.  The  target 

initial  position  is  at  x  =  10000  ft,  y  =  500  ft.  Applied 

lateral  target  acceleration  is  away  from  the  missile. 


6 


in.n,uuuttU  Al  bOVLHNMErn  EXPENSE 


Figure  II-2  Position  Plot  for  Head-on 
Aspect  No  Target  Turn 
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Figure  II-4  Line  of  Sight  Angle  Plot  for  Head-on 
Aspect  No  Target  Turn 
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Fiaure  II-7  Missile  Heading  Plot  for  Head-on 
g  Aspect  Constant  3G  Target  Turn 
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Figure  11-11  Missile  Headings  Plot  for  Head-on 
Aspect  Constant  6G  Target  Turn 
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Figure  11-14  Position  Plot  for  Tail 
Aspect  Constant  No  Target  Turn 
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Figure  11-15  Missile  Heading  Plot  Tail 
Aspect  Constant  No  Target  Turn 
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Figure  11-16  Line  of  Sight  Angle  Plot  Tail 
Aspect  Constant  No  Target  Turn 
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Pigure  11-17  Line  of  Sight  Rate  Plot  Tail 
Aspect  Constant  No  Target  Turn 
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Figure  11-18  Position  Plot  Tail 
Aspect  Constant  3G  Target  Turn 
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Figure  11-19  Missile  Heading  Plot  Tail 
Aspect  Constant  3G  Target  Turn 
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Figure  11-21  Line  of  Sight  Rate  Plot  Tail 
Aspect  Constant  3G  Target  Turn 
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Figure  11-23  Missile  Heading  Plot  Tail 
Aspect  Constant  6G  Target  Turn 
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Figure  13-25  Line  of  Sight  Rate  Plot  Tail 
Aspect  Constant  6G  Target  Turn 
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Figure  IT-27  Missile  Heading  Plot 
Aspect  Constant  No  Target  Turn 
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Figure  11-28 
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Figure  11-30  Position  Plot  Beam 
Aspect  Constant  3G  Target  Turn 
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Figure  11-31  Missile  Heading 
Aspect  Constant  3G  Target 
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Figure  11-33  Line  of  Sight  Rate  Plot  Beam 
Aspect  Constant  3G  Target  Turn 
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Figure  11-35  Missile  Headino  Plot  Beam 
Aspect  Constant  6G  Target  Turn 
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Figure  11-36  Line  of  Sight  Angle  Plot  Beair 
Aspect  Constant  6G  Target  Turn 
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By  applying  the  acceleration  away  from  the  missile  the 
pilot  will  lose  sight  of  the  missile.  This  is  undesirable, 
but  a  turn  into  the  missile  will  help  the  missile,  in  the 
early  stages,  more  than  a  turn  away.  With  a  turn  into  the 
missile,  the  pilot  will  also  lose  sight  of  the  missile, 
during  a  constant  acceleration  turn. 

1 .  OG  Target  Acceleration 

From  the  position  plot,  Figure  II-2,  it  is  seen  that 
proportional  navigation  guidance  is  essentially  the  same  as 
the  direct  path  guidance.  Any  errors  are  due  to 
initialization  of  the  heading  for  the  proportional 
navigation  guidance.  The  pure  pursuit  guidance  missile  and 
the  lead  pursuit  guidance  missile  fly  curvilinear  paths  to 
targetintercept .  The  curve  for  the  lead  pursuit  guidance 
missile  is  less  than  the  pure  pursuit  missile  due  to  target 
lead. 

The  missile  headings  graph,  Figure  II-3,  shows  the 
relative  heading  changes  involved  for  each  missile  guidance. 
After  the  initialization  errors  have  been  corrected,  the 
proportional  navigation  guidance  missile  parallels  the 
direct  path  missile.  The  heading  changes  for  the  lead 
pursuit  are  less  than  for  the  pure  pursuit  guidance  method. 
Large  increases  in  missile  headings  at  the  end  of  the 
intercept  implies  large  lateral  accelerations  are  required 
for  the  missile  to  complete  the  intercept. 

The  line  of  sight  graph.  Figure  II-4,  shows  what 
would  be  expected  for  this  case.  The  proportional 
navigation  guidance  and  direct  path  missile  maintain 
constant  line  of  sight,  approximately,  while  the  line  of 
sight  increases  for  lead  pursuit  and  pure  pursuit  guidance 
methods.  The  large  change  in  line  of  sight  at  the  end  of 
the  intercept  also  correlates  to  a  high  lateral  acceleration 
required  by  the  missile. 

The  line  of  sight  rate  graph,  Figure  11-5,  gives 
some  insight  to  the  control  inputs  to  the  missile  guidance 
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I  subsystem.  The  values  from.  Figure  II-5,  are  the  slew  rated 

&  for  the  sensor  subsystem.  A  positive  slew  rate  is  seen  for 

^  the  prop  nav  guidance  only  at  the  final  stages  of  the 

i  intercept.  Lead  pursuit  and  pure  pursuit  guidance  methods 

j  have  accelerating  positive  slew  rates  throughout  the 

|  intercept.  The  direct  path  missile  has  a  negative  slew 

i  rate,  caused  by  the  missile  speed  advantage  (2500  :  667 

i  ft/sec) . 

S  2.  3G  Target  Acceleration 

I  The  position  graph,  Figure  II-6,  shows  a  curvilinear 

path  for  all  three  guidance  methods.  The  curvature  of  the 
target  flight  path  is  misleading,  because  of  the  axis 
scaling.  The  target  is  maintaining  a  constant  acceleration. 
|  All  three  guidance  methods  appear  to  end  up  in  a  tail  chase. 

The  scaling  is  misleading  again.  Target  heading  change  is 
approximately  50  degrees.  The  proportional  navigation 
missile  impacts  in  the  beam  while  lead  pursuit  will  be  rear 
I  quarter  and  pure  pursuit  will  be  a  tail  chase. 

The  missile  headings  graph.  Figure  II-7,  shows  the 
proportional  navigation  guidance  has  the  lowest  heading 
slope  and  is  approximately  linear  at  the  end  of  the 
intercept.  Pure  pursuit  and  lead  pursuit  guidance  methods 
have  accelerating  missile  heading  slopes  requiring  higher 
missile  acceleration. 

The  line  of  sight  graph.  Figure  II-8,  is  similar  to 
the  Figure  II-4,  proportional  navigation  guidance  method, 
which  has  low  line  of  sight  angles,  slightly  increasing  due 
to  target  acceleration.  Lead  pursuit  and  pure  pursuit 
guidance  methods  have  line  of  sight  angles  which  increase  at 
an  accelerated  rate  throughout  the  intercept.  The  direct 
missile  has  decreasing  line  of  sight.  The  large  negative 
LOS  for  the  direct  missile  at  the  end  of  the  intercepts  is 
caused  by  the  miss  distance  and  heading  initialization. 

Line  of  sight  rates,  Figure  II-9,  correlate  with  the 
line  of  sight  plot,  Figure  II-8.  Line  of  sight  rates  are 
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small  but  show  an  acceleration  at  the  end  of  the  intercept, 
due  to  decreasing  range.  Proportional  navigation  guidance 
methods  are  reducing  the  line  of  sight  angle  while  lead 
pursuit  and  pure  pursuit  increase  the  line  of  sight  angle. 

3 .  6G  Target  Acceleration 

Comparing  the  position  plot,  Pigure  11-10,  with  that 
of  the  3G  case.  Figure  11-6,  similar  statements  can  be  made 
about  all  of  the  missile  paths.  Scaling  is  slightly 
deceiving;  the  target  has  made  approximately  100°  heading 
change.  All  flight  paths  are  curvilinear  with  the 
proportional  navigation  guidance  method  being  the  shorter  of 
the  three  methods. 

Figure  11-11,  shows  the  heading  changes  for  the 
missiles  and  the  smaller  missile  maneuvering  required  for 
the  proportional  navigation  missile.  Figure  11-12  and 
Figure  11-13  show  larger  magnitudes  for  line  of  sight  angle 
and  line  of  sight  rate  than  the  3G  case,  but  follow  the  same 
trends.  The  direct  path  missile  shows  a  reversal  in  line  of 
sight  rate  as  the  target  heading  change  is  greater  than  90°. 

B.  TAIL  ASPECT 


Figure  11-14  through  Figure  11-25  are  the  results  of 
missile  guidance  comparisons  for  tail  aspect  initial 
condition  with  0G,  3G  and  6G  constant  target  acceleration. 
The  missile  begins  at  the  origin  of  the  graph.  The  target 
initial  position  is  X*10000  ft.  and  Y-1000  ft.,  with  an 
initial  heading  of  090,  parallel  to  the  X  axis.  Applied 
target  acceleration  is  directed  into  the  missile, 
perpendicular  to  the  target  heading. 

1 .  0G  Target  Acceleration 

The  position  plot.  Figure  11-14,  shows  the 
proportional  navigation  guidance  missile  flies  a  similar 
path  as  the  direct  path  missile.  The  difference  in  the 
flight  paths  is  due  to  errors  in  initialization.  The 


missile  heading  plot,  Figure  11-15,  shows  that  the  heading 
for  proportional  navigation  guidance  and  direct  path 
missiles  are  parallel  after  the  initialization  errors  are 
corrected.  Pure  pursuit  and  lead  pursuit  guidance  methods 
have  continually  changing  headings  with  accelerating  slopes 
at  the  final  stage  of  the  intercept. 

Line  of  sight  angles  for  the  proportional  navigation 
guidance  and  direct  path  are  approximately  constant  and 
equal  to  the  initial  line  of  sight  angle,  giving  a  constant 
bearing  decreasing  range  trajectory  as  seen  by  the  target. 
The  line  of  sight  angle  for  pure  pursuit  and  lead  pursuit 
guidance  decrease,  but  non  linearly,  as  seen  in  Figure  11-16 
and  Figure  11-17 . 

2 .  3G  Target  Acceleration 

With  target  acceleration,  all  three  missiles  fly  a 
curvilinear  path.  The  proportional  navigation  guidance 
method  has  the  shortest  flight  path,  as  seen  in  Figure  Il¬ 
ls.  Proportional  navigation  guidance  gives  a  linear  heading 
change,  as  seen  in  Figure  11-19.  Pure  pursuit  and  lead 
pursuit  guidance  methods  give  higher  heading  slopes  when  the 
target  applies  lateral  acceleration  as  compared  to  the  OG 
heading  plot,  Figure  11-15. 

Line  of  sight  angle  changes  are  small  for  propor¬ 
tional  navigation  guidance,  as  shown  in  Figure  11-20.  Pure 
pursuit  and  lead  pursuit  guidance  have  decreasing  line  of 
sight  angles,  with  corresponding  decreasing  line  of  sight 
rates,  as  seen  in  Figure  11-20  and  Figure  11-21.  Line  of 
sight  rates  increase  for  proportional  navigation  guidance, 
as  would  be  expected  from  the  path  the  missile  flies.  The 
direct  path  gives  both  a  nonlinear  line  of  sight  angle  and 
line  of  sight  rate  throughout  the  intercept. 

3 .  6G  Target  Acceleration 


When  the  target  acceleration  is  increased,  flight 
paths  have  a  larger  curvature,  as  seen  in  Figure  11-22.  The 
scaling  gives  some  distortion,  the  target  has  gone  through 
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approximately  100°  of  heading  change.  The  missile  heading 
changes,  as  per  Figure  11-23,  are  similar  to  those  observed 
for  the  6G  head-on  aspect,  Figure  11-11. 

The  line  of  sight  angle  and  line  of  sight  rate  are 
larger  for  an  increase  in  lateral  target  acceleration,  as 
seen  by  comparing  Figure  11-24  and  Figure  11-25  with  the  3G 
case,  Figure  11-20  and  Figure  11-21.  The  line  of  sight  and 
line  of  sight  rate  for  the  direct  path  have  a  very  large 
slope  at  the  final  intercept  due  to  effects  of  decreased 
range.  The  reversal  of  line  of  sight  rate  for  the  direct 
path  in  Figure  11-25  is  where  the  target  heading  change  is 
90°. 


C.  BEAM  ASPECT 

Figures  11-26  through  11-37  are  the  result  of  missile 
guidance  comparisons  for  beam  aspect  initial  conditions  with 
0G,  3G  and  6G  constant  target  accelerations.  The  missile 
begins  at  the  origin  of  the  graph.  The  target  initial 
position  is  X«15000,  Y=0.  Applied  acceleration  is  directed 
into  the  missile. 

1 .  0G  Target  Acceleration 

As  in  the  two  previous  cases,  with  no  lateral  target 
acceleration,  proportional  navigation  guidance  and  direct 
path  missiles  have  similar  flight  paths.  Pure  pursuit  and 
lead  pursuit  guidance  have  curvilinear  flight  paths,  as  seen 
in  Figure  11-26.  Heading  changes  are  small  for  proportional 
navigation  guidance  and  direct  flight  path  missiles  but  not 
zero  as  what  might  be  inferred  from  Figure  11-27,  because  of 
scaling.  The  heading  change  for  pure  pursuit  and  lead 
pursuit  guidance  is  accelerating  throughout  the  flight  time 
with  the  intercept  ending  in  a  tail  chase. 

Line  of  sight  and  line  of  sight  rate,  Figure  11-28 
and  Figure  11-29,  are  similar  to  the  two  previous  cases,  for 
no  target  acceleration  and  the  analysis  is  the  same. 


47 


2 .  3G  Target  Acceleration 

For  the  beam  aspect  initial  condition,  when  lateral 
acceleration  is  applied,  the  difference  between  flight  paths 
for  proportional  navigation  guidance  and  direct  path  is 
opposite  from  the  two  previous  cases  for  lateral 
acceleration.  The  proportional  navigation  guidance  missile 
flight  path  is  on  the  opposite  side  of  the  direct  path  from 
pure  pursuit  and  lead  pursuit  guidance  flight  paths,  Figure 
11-30,  with  the  opposite  curvature.  Headings  for 
proportional  navigation  guidance  continually  decrease  while 
pure  pursuit  and  lead  pursuit  guidance  increase. 

Differences  between  the  methods  are  enhanced  by  the 
line  of  sight  angle  and  line  of  sight  rate  plots  in  Figure 
11-32  and  Figure  11-33.  Proportional  navigation  guidance 
decreases  line  of  sight  while  the  others  have  an  increasing 
line  of  sight  and  appropriate  line  of  sight  rate. 

3 .  6G  Target  Acceleration 

Increasing  the  target  lateral  acceleration  magnifies 
the  flight  path  differences  between  the  guidance  methods. 
As  has  been  seen  from  the  previous  cases,  the  larger  lateral 
acceleration  increases  the  magnitudes  of  the  values  for 
Figure  11-33  through  Figure  11-37,  compared  with  similar 
graphs  from  the  other  cases,  but  the  trends  remain  the  same. 
Proportional  navigation  guidance  parameters  have  smaller 
changes  than  pure  pursuit  and  lead  pursuit  guidance,  with 
parameters  generally  decreasing  instead  of  increasing  for 
pure  pursuit  and  lead  pursuit  guidance. 

D.  CONCLUSIONS 

For  scenarios  with  no  applied  target  lateral 
acceleration  the  proportional  navigation  missile  is  the  same 
as  the  direct  missile.  The  pure  pursuit  and  lead  pursuit 
missiles  finish  in  a  'tail  chase'  where  a  missile  speed 
advantage  is  required  to  complete  the  intercept.  The  line 
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of  sight  remains  constant  for  the  proportional  navigation 
missile  but  increases  with  the  pursuit  missiles.  The  line 
of  sight  rate  increases  with  decreased  range  for  the  pursuit 
missiles  but  is  zero  for  the  proportional  navigation  and 
direct  missiles. 

When  target  lateral  acceleration  is  applied,  there  is  a 
deviation  between  the  direct  and  the  proportional  navigation 
missiles.  Since  the  target  is  turning  into  the  missile,  the 
line  of  sight  angle  decreases  at  an  accelerated  rate  as 
range  decreases.  The  proportional  navigation  missile 
accounts  for  the  change  of  line  of  sight  by  turning  into  the 
target.  The  pursuit  missiles  fly  a  tail  chase  profile  with 
higher  line  of  sight  accelerations  due  to  the  target  turn. 

For  the  direct  missile,  when  target  acceleration  is 
applied,  the  line  of  sight  is  not  constant  and  the  rate  of 
change  depends  on  the  applied  acceleration.  Implementation 
of  a  direct  missile  is  impossible  because  the  parameters 
used  to  guide  the  missile  are  dependent  on  the  target  flight 
path. 

If  an  optimum  missile  is  to  be  designed,  proportional 
navigation  guidance  is  the  closest  to  an  "ideal”  missile. 
The  better  the  proportional  navigation  missile  can 
compensate  for  the  effects  of  the  target  acceleration,  the 
closer  to  "ideal"  the  missile  will  become. 
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III.  TARGET  MODEL 


A  complete  target  model  for  use  in  computer  simulation 
is  very  involved,  time  consuming  and  computer  intensive.  To 
simplify  target  simulation  the  target  flight  profile  is 
based  on  the  fact  that  the  missile  sees  only  the  effects  of 
the  target  command  inputs  and  resultant  flight  path.  The 
target  model  was  simplified  to  include  only  the  flight 
profile  desired  and  not  be  concerned  with  the  full  target 
modeling.  A  constant  speed,  constant  acceleration  target  is 
assumed  for  the  simulation.  A  variable  speed,  variable 
acceleration  target  can  be  added  at  a  later  time.  The 
missile  simulation  estimates  and  predicts  target  parameters 
of  range,  range  rate,  range  acceleration,  bearing,  bearing 
rate,  and  bearing  acceleration.  Therefore,  for  proper 
evaluation  of  the  missile  guidance  and  missile  flight 
profiles  the  target  parameters  in  missile  coordinates  for 
all  these  parameters  must  be  computed. 

Complete  analysis  of  target  motion  is  obtained  from  a 
three  dimensional  derivation,  but  insight  can  be  gained  from 
two  dimensional  modeling.  Three  dimensional  flight  profiles 
are  easily  implemented  on  the  computer  but  graphic  display 
of  the  results  are  difficult.  Two  dimensional  displays  are 
easier  to  implement  and  comprehend.  A  two  dimensional 
target  model  is  assumed. 

A  target  can  accelerate  at  values  ranging  from  negative 
maximum  instantaneous  acceleration  to  positive  maximum 
instantaneous  acceleration,  A(max  inst) .  A(max  inst)  is 
defined  as  the  aerodynamic  acceleration  given  by  the  maximum 
deflection  of  control  surfaces.  A (max  inst)  is  dependent  on 
airspeed  and  air  density.  High  speeds  and  low  altitudes 
produce  the  highest  instantaneous  accelerations.  Maximum 
sustained  acceleration,  A(max  sust) ,  is  defined  as  the 
aerodynamic  acceleration  to  maintain  constant  airspeed  and 
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constant  altitude,  at  full  thrust.  A (max  sust)  can  be 
exceeded  but  must  be  compensated  for  by  a  reduction  in 
airspeed  or  altitude. 

In  three  dimensional  maneuvering  cross  coupling  exists 
between  horizontal  and  vertical  angle  and  angle  rate 
components  of  target  velocity  and  target  acceleration. 
Applied  accelerations  and  velocity  changes  in  one  direction 
will  affect  the  parameters,  seen  by  a  missile,  in  the  two 
other  directions. 

Thrust  capabilities  have  a  direct  correlation  to  A (max 
sust)  and  the  airspeed  of  an  aircraft.  An  aircraft  with 
higher  thrust  can  maintain  a  higher  speed  and  compensate  for 
drag  induced  by  the  applied  acceleration.  A  modern  aircraft 
with  a  relatively  high  thrust  to  weight  ratio  will  have  a 
very  high  A(max  sust)  which  is  close  to  A(max  inst) . 

Airspeed  is  a  key  element  for  maneuverability  and 
survivability.  Tactics  incorporate  optimum  techniques  for 
maintaining  airspeed  or  recovering  lost  airspeed.  Pilots 
learn  to  compensate  for  limitations  of  A (max  sust)  by 
intentionally  decreasing  altitude  and  use  the  effects  of 
gravity  to  maintain  airspeed  when  lateral  acceleration  is 
applied.  Another  technique  is  to  apply  the  lateral 
acceleration  required  to  perform  a  maneuver  then  to  reduce 
the  acceleration,  allowing  excess  thrust  to  restore  the 
airspeed  and  altitude  lost  during  the  maneuver.  There  is  a 
recovery  time  for  the  thrust  to  restore  the  lost  energy,  so 
to  aid  in  restoring  airspeed,  a  pilot  will  normally  go  to 
zero  acceleration,  reducing  any  induced  drag,  effectively 
increasing  the  aircraft  thrust.  This  maneuver  causes  a  loss 
in  altitude  due  to  gravity  but  improves  airspeed 
restoration. 

The  probability  of  aircraft  acceleration  is  used  to 
determine  parameters  for  the  target  model  used  in  missile 
designs  and  simulations.  Figure  III-l  shows  a  typical 
acceleration  probability  graph  used  in  missile  design.  The 
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figure  does  not  account  for  pilot  tendencies  nor  the 
difference  between  A (max  sust)  and  A (max  inst) .  The  graph 
assumes  that  a  pilot  will  maneuver  primarily  at  zero 
acceleration ,  straight  and  level ,  or  maximum  acceleration, 
for  a  turn,  with  some  probability  for  any  other  possible 
acceleration.  The  design  engineer  assigns  probabilities  for 
the  impulse  functions  at  zero  acceleration  and  at  maximum 
acceleration  depending  on  the  type  of  target  aircraft.  A 
large  bomber  may  have  an  A (max  sust)  half  that  of  a  fighter 
aircraft,  with  less  probability  of  turning  than  flying 
straight  and  level. 


Figure  III-l  Probability  of  Aircraft  Acceleration  [Ref.  2,3] 

A  proper  target  maneuver  model  should  include  some  pilot 
tendencies  and  known  tactics.  A  pilot  is  not  always  able  to 
move  the  control  surface  to  a  precise  location  to  cause  a 
precise  acceleration  at  an  optimal  time.  A  pilot  will  move 
the  control  surface,  judge  the  acceleration  induced  then 
move  the  control  surface  to  achieve  a  desired  acceleration. 
The  feel  a  pilot  receives  from  the  'stick'  is  a  prime 
feedback  source  to  allow  the  pilot  to  set  the  desired 
acceleration.  The  more  force  the  pilot  applies,  the  more 
the  control  surface  moves,  and  the  higher  is  the 
acceleration.  Modern  aircraft  may  use  computers  to  achieve 
the  commanded  acceleration,  reducing  any  pilot  induced 
errors  on  input. 
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A.  TACTICS  FOR  MISSILE  DEFENSE 


When  a  missile  is  fired  at  the  target  the  battlefield 
scenario  changes  to  an  immediate  survival  situation.  If  the 
missile  is  undetected,  the  acceleration  probabilities  of 
Figure  III-l  may  be  an  adequate  target  model.  If  a  pilot 
sees  the  missile,  pilot  reactions  will  change  the 
probabilities.  How  the  probabilities  change  may  be  of 
consequence  to  the  missile  guidance.  An  optimum  missile 
design  may  be  able  to  use  pilot  tendencies  to  increase 
missile  performance.  The  overall  acceleration  probability 
from  Figure  III-l  is  a  zero  mean  function  with  a  variance, 
o2,  dependent  on  the  probabilities  assigned.  Reference  2 
discusses  obtaining  the  parameters  for  the  target  accel¬ 
eration  probability  model  of  Figure  III-l. 

When  the  pilot  imposes  a  missile  defense,  the  overall 
acceleration  may  not  be  zero  mean,  nor  maintain  the  same 


variance.  To  account  for  the  changes  in  acceleration 
probability  a  function  similar  to  Figure  III-2  might  be 
used.  This  model  accounts  for  some  variance  to  the 
acceleration  which  the  pilot  is  trying  to  achieve  centered 
around  A(max  sust) .  If  the  pilot  is  trying  to  achieve  A(max 
inst) ,  it  is  assumed  he  will  be  decreasing  airspeed  and 
reducing  actual  acceleration  until  the  applied  acceleration 
is  decreased  to  A(max  sust)  or  below.  A  smart  pilot  will 
either  fly  at  a  maximum  acceleration  or  at  zero 
acceleration,  increasing  the  aircraft  maneuverability. 

Last  ditch  maneuvers  are  performed  at  A (max  inst)  to 
avoid  the  missile,  neglecting  any  adverse  effects  of 


applying  acceleration,  in  order  to  increase  survivability. 


If  the  last  ditch  maneuver  is  performed  too  soon. 


acceleration  is  decreased,  due  to  loss  of  airspeed,  negating 


Figure  III-2  Probability  of  Acceleration 

If  pilot  reaction  is  taken  as  Gaussian  when  trying  to 
achieve  a  desired  acceleration,  the  overall  acceleration 
probability  will  be  Gaussian  with  a  non-zero  mean  and  a 
variance  dependent  on  the  combination  of  the  two  Gaussian 
terms.  The  probability  assigned  to  each  term  determines  the 
mean  and  variance. 

Missile  defense  includes  placing  the  missile  on  the 
beam,  to  utilize  the  largest  acceleration  vector,  with 
maneuvers  made  out  of  phase,  out  of  plane  with  the  missile. 
The  largest  acceleration  component  comes  from  the  elevator, 
perpendicular  to  the  wings.  By  placing  the  wings  parallel 
with  the  plane  of  the  missile,  the  largest  acceleration 
component  is  used  to  create  the  largest  missile  corrections, 
perpendicular  to  the  plane  of  the  missile.  The  plane  of  the 
missile  is  defined  by  three  points:  target  position, 
missile  position  and  the  projected  impact  point. 

A  graph  of  a  target  acceleration,  while  performing 
missile  evasion,  might  look  like  Figure  II1-3.  The  pilot 
commands  A(max  sust)  or  A(max  inst)  for  a  short  time,  then 
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reduces  the  acceleration  to  zero  to  regain  lost  energy, 
before  applying  the  acceleration  again.  This  process  may  be 
repeated  2,3  or  more  times  during  the  missile  flight  time. 
The  graph  attempts  to  incorporate  transient  response  induced 
by  system  time  delays,  transient  responses  of  the  control 
surfaces  and  pilot  tendencies  for  maneuvering  and  control  of 
the  target. 


TIME  (sec)' 


Figure  III-3  Target  Acceleration 

The  resultant  average  acceleration  is  non-zero,  with  a 
non-zero  variance.  The  mean  and  variance  of  Figure  III-3 
can  be  estimated  by  the  parameters  assumed  in  Figure  III-2. 

Figure  III-l  and  Figure  III-2  show  total  aircraft 
acceleration.  The  parameters  as  seen  by  the  missile,  in 
antenna  coordinates,  will  vary  depending  on  the  three 
dimensional  maneuver  employed.  The  missile  tracking 

subsystem  must  be  designed  to  handle  the  maximum 
acceleration  possible  for  each  of  the  orthogonal  components 
of  the  reference  frame. 

Proportional  navigation  missiles  compensate  for  the 
applied  target  acceleration  by  decreasing  the  line  of  sight 
rate  induced  by  the  change  of  target  velocities.  With  a 
good  target  model  and  detection  of  maneuvering  effects,  the 
missile  guidance  can  predict  target  motion  and  position. 


IV. 


SYSTEM  MODEL 


From  the  section  on  the  ideal  missile,  it  was 
ascertained  that  target  parameters  are  not  constant  if 
lateral  acceleration  is  applied.  This  thesis  will  attempt 
to  incorporate  as  much  sensor  information  as  is  available  in 
defining  the  system  models  for  an  optimum  missile.  Current 
radars  allow  the  measurement  of  range,  range  rate  and  off 
boresight  bearing  error.  Measurements  taken  by  the  radar 
are  referenced  to  a  radar  axis  system.  In  order  for  the 
missile  to  use  the  information  supplied  by  the  radar,  a 
common  reference  frame  must  be  established. 

A.  COORDINATE  SYSTEM 

Each  entity  in  the  missile-target  intercept  problem  has 
its  own  coordinate  system.  The  overall  geometry  as  seen 
from  an  "eye  in  the  sky"  would  view  it  in  space  coordinates. 
An  observer  on  the  ground  would  view  it  in  earth 
coordinates.  The  launching  aircraft,  missile  and  target 
aircraft  will  view  it  in  an  individual  coordinate  system, 
referenced  to  that  specific  platform.  Trying  to  equate  each 
coordinate  system  is  not  an  easy  task  but  one  which  is  done. 

By  use  of  Euler  angles  any  reference  coordinate  system  can 
be  related  with  Earth  coordinates.  By  use  of  a  directional 
cosine  matrix  transformation  any  reference  coordinate  system 
can  be  transformed  to  another  reference  coordinate  system 
[Ref .  3] . 

The  missile  is  concerned  with  flying  to  a  point  in  space 
that  will  hopefully  be  occupied  by  the  target  at  the 

completion  of  the  intercept.  The  object  is  to  guide  the 

missile  to  the  proper  point  where  the  target  will  be.  The 
missile  is  concerned  with  its  coordinate  system  and  not 

that  of  the  target.  But  on  the  missile  itself  there  are 
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various  coordinate  system  reference  points.  Each  sensor  has 
its  own  location  on  the  missile  and  where  it  is  mounted  is 

its  reference  point.  Any  moving  sensor,  like  the  antenna, 

will  have  its  special  reference  coordinate  system.  Missile 
parameters  are  normally  referenced  to  the  missile  body  frame 
of  reference  while  target  parameters  are  referenced  to  the 
antenna  frame  of  reference.  While  very  complicated,  the 
frames  of  references  can  be  transformed  and  equated. [Ref .  3] 
To  simplify  simulation  and  evaluation  of  desired 
parameters,  an  inertial  frame  of  reference  will  be  used 
which  is  centered  at  the  radar  antenna  location.  This 
simplification  will  aid  in  better  evaluation  of  the  effects 
of  the  target  parameters  and  the  missile  guidance  without 
encumbrance  of  transformation  errors.  Although  the 

simplification  assumes  ideal  missile  parameters,  time  delays 
can  be  incorporated  later  to  account  for  first  order 

modeling  of  the  missile. 


r 


B. 


EQUATIONS  OF  MOTION 


L 


In  cartesian  coordinates  missile  and  target  motion  is 
described  by  the  standard  motion  equation: 


X  (T)  «  X0  + 


‘T 

,  o 


X(t)  dt  + 


‘  PT 

.  Jo 


X(t)  dt  dt 


(7) 


The  equation  is  based  on  a  fixed  reference  point.  The 
orthogonal  directions  (Y  and  Z)  will  have  the  same  equation. 

The  antenna  frame  of  reference  uses  polar  coordinates 
which  have  the  equations: 


r(t)  *  r*  +  r(t)  dt  + 


I/(t>  dt  +  ||Iri 

e(t)  dt  +  jj^  ©( 

II iU)  dt  +  ||I  4 


(t)  dt  dt 


e(t)  »  e»  +  |  8 ( t )  dt  + 

4>(t)  =$©  +  I  $(t)  dt  + 


(t)  dt  dt 
4»"(t)  dt  dt 


(8) 

(9) 

(10) 
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where  r(t)  =  the  radial  component  of  motion 
©(t)  *  horizontal  component  of  motion 
♦  (t)  *=  vertical  component  of  motion 

When  the  reference  point  is  not  fixed,  extra  terms  and 
cross  coupling  are  introduced  into  equation  dynamics.  The 
coriolis  equation  accounts  for  the  moving  reference  point. 
Depicted  in  Figure  IV-1  a  change  in  the  vector  R  is 
accounted  for  by  both  changes  in  the  magnitude  and  the 
rotation  effects  by  the  moving  reference  point.  Using  the 
terms  as  defined  in  Figure  IV-1  we  can  obtain  the  necessary 
equations  to  find  r(t),  d(t)  and  4>(t).  For  simplicity,  only 
the  derivation  for  r  and  ©  are  shown  with  r  and  4> 
relationships  being  a  duality  of  derivation  of  r  and  ©.  The 
simulation  of  Section  VI  will  be  two  dimensional. 


Figure  IV-1  Rotating  Vector  Diagram 
The  coriolis  equation  to  relate  the  time  rate  of  change 
of  the  R  vector  to  the  rate  of  change  in  the  r  direction  and 
the  angular  rate  of  motion  is: 

R  »  r  +  wxR  (11) 

where  R  =  the  directional  vector 


r  =  the  magnitude  of  the  directional  vector 
w  =  the  angular  rate  of  motion 


The  total  change  of  the  vector  R  is  the  sum  of  the  change  in 
the  magnitude  of  R  due  to  changes  along  the  original  vector 
R  given  by  r  and  the  angular  rotation  due  to  the  moving 
coordinate  frame,  given  by  wxR. 

Utilizing  the  general  rule  for  differentiating  a  vector1 
an  expression  is  obtained  for  the  acceleration  of  the  R 
vector. 


R  « 

r 

+ 

wxr 

+ 

wxr 

+  wxR 

(12) 

«  • 

R  * 

•  ■ 

r 

+ 

wxr 

+ 

wxr 

+  wxr  +  wxwxr 

(13) 

•  • 

R  = 

r 

+ 

wxr 

+ 

2 (wxr)  +  wxwxr 

(14) 

where  R  *  directional  vector 

r  =  the  acceleration  of  the  magnitude  of  R 
wxr  and  wxr  are  perpendicular  to  the  R  vector 
wxwxr  is  centrifugal  acceleration 


This  equation  gives  the  cross  correlation  of  range  and 
angle  to  implement  in  a  second  order  model.  Applying  the 
rule  of  differentiating  a  vector  again  will  yield  the 
equations  for  a  third  order  model. 


^ 

+  wxwxr 


+  w  xr  +  wxr  +  2 (wxr) 
+  wxwxr  +  wxR 


+  2 (wxr  )  +  wxwxr 

(15) 


d (R  )  *  r*  +  wxr  +  wxr  +  2 (wxr)  + 
d  t  .  .  ... 

+  wxwxr  +  wxwxr  +  wxr  +  wxwxr  + 

R  *  r  +  3  (wxr)  +  3  (wxr  )  +  w  xr 
+  2 (wxwxr)  +  3 (wxwxr)  +  wxwxwxr 


2 (wxr  )  +  wxwxr 

wx2(wxr)  +  wxwxwxr 

(16) 


+  wxwxr 


(17) 


•  •  • 

where  R  is  the  acceleration  jerk  of  vector  R 


1  rule  for  differentiating  a  vector 
^  +  wxA 

the  total  derivative  is  the  sum  of  the  time  rate  of  change 
of  the  vector  and  the  rotation  of  the  vector. 
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r  is  the  magnitude  of  R 

wxr  is  perpendicular  to  R 

wxwxr  is  in  the  negative  direction  of  R 

wxwxwxr  is  perpendicular  to  R 

C.  SECOND  ORDER  MODEL 


Beginning  with  equations  8-10, 
space  system  can  be  set  up  which  would 
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a  second  order  state 
have  the  form: 


(18) 

(19) 


(20) 


The  range  portion  of  the  second  order  system  may  be 
accomplished  totally  by  the  radar,  since  no  other  subsystems 
require  the  information.  The  radar  receiver  is  designed  to 
track  the  target  in  the  radial  direction  without  the  need 
for  an  additional  filter. 

The  more  difficult  state  equations  to  implement  in  the 
missile  are  the  angular  directions.  The  second  order,  time 
invariant,  constant  velocity,  zero  acceleration,  state 
feedback  model  makes  the  tracking  much  easier.  The 
continuous  system  model  can  be  given  in  a  time  derivative 
form  as : 


x’  +  ki  -X  +  k*  X  «  0 


(21) 


If  the  term  given  in  the  equation  as  X  is  actually  the  error 
of  the  angular  position  then 

X  =  (line  of  sight  angle  -  antenna  position  angle) 
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where  X  can  be  directly  measured  by  the  antenna.  The  values 
of  the  k's  in  the  time  derivative  equation  depend  on  the 
designer  and  the  response  desired.  For  the  simulations  of 
the  proportional  navigation  missiles  of  Section  II  and 
Section  VI,  ki  =  20  and  k2  =  100  were  used.  These  constants 
give  a  response  time  constant  of  0.1  sec. 

The  use  of  the  coriolis  equation  to  derive  the  second 
order  model  gives  a  time  varying  solution.  Implementing 
time  varying  equations  are  difficult  and  often  avoided  by 
using  a  time  invariant  system  and  state  feedback  to  cancel 
errors.  The  time  variant  space  state  model  derived  from 
equation  14  is  shown  below: 

R  *  r  +  wxr  +  2(wxr)  +  wxwxr  (22) 


Separating  into  orthogonal  components  of  radial  and 
transverse  with  scalar  multiplication: 


a*  =  r  +  wwr  (23) 

ai  =  wr  +  2(wr)  (24) 

Rearranging  into  equations  to  implement  into  a  system: 


r  *  =  -© 2r  +  as  (25) 

6  *  =  -2©r  +  aT  (26) 

r  r 

The  state  space  model  is  difficult  to  represent  unless 
divided  into  two  channels  with  cross  coupling  given  in  the  A 
matrix. 
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D.  THIRD  ORDER  MODEL 

The  time  invariant  model,  derived  from  equations  similar 
to  those  deriving  the  time  invariant  second  order  model 
(equations  18  and  19) ,  in  state  space  form  is  given  by: 
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Ranging  may  be  accomplished  by  the  radar  receiver,  as  in 
the  second  order  model,  for  similar  reasons.  Tracking 
angular  positioning  requires  knowledge  of  the  target  angular 
acceleration  as  well  as  angular  velocity.  A  filter  is 
normally  used  to  maintain  a  track  of  the  target  angular 
parameters.  Common  filters  are  a-p ,  Weiner,  and  Kalman.  As 
compared  in  Reference  4,  the  Kalman  is  the  best  filter 
suited  for  air  to  air  missiles,  but  also  the  most  costly  to 
implement.  For  the  time  invariant  third  order  model  a 
simple  constant  gain  Kalman  Filter  can  be  used.  The  Kalman 
Filter  will  be  discussed  in  the  next  section. 

The  tine  variant  third  order  model  is  obtained  from  the 
second  derivative  of  the  coriolis  equation  derived  in  the 
previous  section. 

R  =  r  +  3(wxr)  +  3(wxr)  +  w  xr  +  wxwxr 

+  2 (wxwxr)  +  3 (wxwxr)  +  wxwxwxr  (30) 

The  R  term  is  the  change  in  the  acceleration  of  the 
vector  or  a  "jerk"  term,  a  simple  comprehension  is  the  rate 
at  which  the  pilot  applies  the  commanded  acceleration. 
Separating  the  equation  into  radial  and  tangential  terms, 
the  two  orthogonal  scalar  equations  are: 

a»  =  r  3w*r  -  3wwr  (31) 
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kt  -  3wr  +  3wr*  +  w  r  +  w3  r  (32) 

Converting  the  equation  into  primary  coordinate  axis, 
the  equations  obtained  are: 

*r*«  at  +  30 *r  +  3©(0*)r  (33) 

'©*-  at  -  3 (0 ‘ ) r  -  30 (r ‘ )  -  O3  (34) 

r  r 


The  disassociated  space  state  model  looks  like: 
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(36) 


It  is  readily  observed  that  a  very  high  cross  coupling 
of  the  radial  and  transverse  components  exists.  The  range, 
range  rate  and  range  acceleration  are  required  to  adequately 
compute  the  angular  acceleration.  The  angular  velocity  and 
acceleration  is  required  in  computing  the  range 
acceleration.  All  of  these  quantities  are  time  varying 
requiring  a  time  varying  filter  to  implement  this  model. 

The  cubic  term  of  angle  rate  in  the  angle  channel  is 
insignificant  compared  to  the  other  terms  and  is  neglected. 
The  second  order  model  uses  the  simplifying  assumption  of 
constant  velocity  and  constant  acceleration.  For  the  full 
third  order  model,  no  simplifying  assumptions  will  be  made. 
This  third  order  model  should  account  for  all  of  the  cross 
coupling  between  the  bearing  and  angle  channels. 

A  Kalman  Filter  can  be  employed  to  track  the  target  in 
both  range  and  bearing  to  implement  the  time  variant  third 
order  model.  The  Kalman  Filter  will  be  discussed  in  the 
next  section. 
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V. 


KALMAN  FILTER 


Given  a  system  model,  where  the  plant  can  be  modeled  by 
a  set  of  first  order  differential  equations  and  the  output 
can  be  measured,  a  set  of  state  equations  can  be  defined 
similar  to: 

X  =  AX  +  BU  +  W  (37) 

Y  «  C  X  +  V  (38) 

where  X  is  the  state  vector 

Y  is  the  system  output  vector 

A , B , C  and  D  are  matrices 

U  is  the  system  input 

W  is  plant  disturbances 

V  is  measurement  noise 

The  system  can  be  modeled  in  discrete  time  as: 

X (k+1 )  *  4>  X  (k)  +  T  U(k)  +  W  (k)  (39) 

Y (k)  *  H  X (k)  +  V(k)  (40) 

A  Kalman  Filter  is  the  best  filter  to  track  the  output 
of  a  discrete  system  [Ref.  3].  The  Kalman  Filter  equations 
are  given  as: 

X(“|k)  =  X(k|k-i)  +  G(k)*fLY(k)  -  Y  (k  |  k  -  i  )  "|  (41) 


X(k4  1  |k  )  =  $-X(k  |k  )  +  T-U(k)  (42) 

Y(k4  1  |k  )  =  H •  X (k  *  1  |k  )  (43) 

A 

where  X(k  |k )  =  the  state  estimate  at  time  k  given 

information  through  time  k. 

X(k41  |k)  =  the  state  estimate  at  time  k+1  given 
information  through  time  k. 

Y(k41  |k )  =  the  output  estimate  at  time  k+1 
given  information  through  time  k. 

=  the  filter  gain  at  time  k. 


G(k) 


For  linear,  time  invariant  systems,  the  ♦  and  r  matrices 

are  easy  to  calculate  and  follow  directly  from  the  state 

At  f  At 

space  model,  where  ♦  *  e  and  r  *  e  dt.  For  non  linear 
systems,  an  extended  Kalman  filter  can  be  used.  For  the 
extended  Kalman  filter,  the  $  and  r  matrices  are  linearized 
about  the  projected  operating  point.  One  method  of 

estimating  the  linearization  is  to  take  the  partial 
derivatives  of  the  non  linear  state  space  matrices: 

(44) 

(45) 
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The  gain  matrix  G(k)  will  vary  with  the  parameters  of 
the  filter.  G(k)  is  the  weighting  factor  for  the  system 
error.  The  solution  to  the  filter  gain  G(k)  requires  the 
solution  of  Riccati  equations: 

G  (k)  =  PC  |  k  -  i  )  HT  [  H  PC  |k-i  >  H*  +  R(k)]  1  (46) 

PC  |k-i  )  «  4>  PC  |k-i  )  *T  +  6  Q  6T  (47) 

PC|k)  =  PC|k-i)  -  G  (k)  *  PC|k-i)  (48) 

where  G(k)  =  Kalman  Filter  gain  at  time  k 

PC|k-i)  =  Covariance  of  predicted  estimate 
R(k)  =  measurement  covariance  matrix,  EfWH 
Q(k)  =  maneuver  covariance  matrix,  ElUUM 
PC|k)  =  Covariance  of  filtered  estimate 
5  =  maneuvering  weighting  matrix 

A  constant  gain  matrix  can  also  be  used  in  the  Kalman 
filter.  Instead  solving  the  full  Riccati  equations  for  each 
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change  of  variables,  a  constant  value  is  used  throughout  the 
problem.  A  constant  gain  matrix  will  simplify 
implementation  of  the  Kalman  filter. 

One  implementation  of  the  third  order  model,  as 
discussed  in  the  previous  section,  is  to  model  the  system  as 
linear,  time  invariant,  given  by  the  space  state  model: 


•  u 

i.  ... 

•  • 

r 

•  ■  • 

r 
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e 

e‘ 

e’ 
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0 

0 
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0 
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0 

0 

0 
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0  0  0  0 

r 

10  0  0 

r 

0  0  0  0 

r 

0  0  10 

e 

0  0  0  1 

e 

0  0  0  0 

e’ 

(49) 


are : 


The  Kalman  Filter  equations  for  this  third  order  model 

A  A  A 

X(k  |k  )  =  X(k  |k-i  )  +  G(k  )  •  [  Y(k)  -  Y(“  |k-i  )  ]  (50) 

A  A 

X(k4  »  |k  )  =  4>*X  (k  |k  )  (51) 


Y(ki  1  |k  )  =  H *X  (k  4  1  |k  ) 


(52) 


where  X  = 


r 

r 

r 

e 

e 

e 


and 


H  = 


1  0  0  0  0  0 
0  1  0  0  0  0 
0  0  1  0  0  0 
0  0  0  1  0  0 
0  0  0  0  1  0 
0  0  0  0  0  1 


Using  a  Kalman  Filter  on  this  third  order  model  is  very 
simple  and  requires  few  on  line  calculations.  The  gain 
matrix  can  be  considered  either  constant  or  time  varying. 
If  time  varying  gains  are  used,  they  can  be  computed  off¬ 
line  and  stored  in  memory.  The  Filter  then  utilizes  the 
precomputed  gain  schedule  and  can  select  a  gain  depending  on 
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the  accuracy  of  the  filter  at  that  time.  If  a  maneuver 
causes  the  filter  to  loose  accuracy  then  a  higher  gain  term 
can  be  utilized.  If  constant  gains  are  used  then  they  must 
be  high  enough  to  compensate  for  any  maneuver  the  target 
might  make.  A  high  gain  matrix  will  make  the  missile  more 
responsive  to  any  unwanted  noise  terms  in  the  system  since 
the  missile  cannot  distinguish  a  noise  input  from  a  target 
maneuver . 

As  discussed  previously,  a  Kalman  Filter  is  not  required 
for  the  range  channel.  The  radar  can  measure  range  and 
range  rate  directly.  Since  the  actual  values  of  the  range 
channel  are  not  used  by  any  other  elements  of  the  guidance 
subsystem,  the  radar  is  able  to  maintain  its  own  tracking  of 
the  target  in  the  center  of  the  range  gate,  which  has  no 
consequences  on  the  rest  of  the  missile  guidance.  Some 
noise  information  can  be  gained  when  estimating  the  range 
channel  with  a  Kalman  Filter.  A  Kalman  Filter  is  used  for 
the  range  channel  in  the  simulation  of  Section  VI  for 
completeness  of  simulation  and  practical  experience. 

A  second  implementation  of  the  third  order  model  is  by 
using  the  equations  obtained  through  the  coriolis  equations. 
The  disassociated  state  space  model  is  given  as: 
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(53) 


(54) 


Two  possible  ways  to  implement  the  Kalman  filter  are  to 
create  an  extended  Kalman  filter  by  linearizing  the  <t> 
matrix,  reducing  the  time  dependence  and  cross  coupling  of 
the  range  and  bearing  channel,  or  keep  the  cross  coupling 
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components  and  have  the  Kalman  Filter  maintain  the  values  of 
the  time  varying  ♦ .  If  ♦  is  linearized,  some  of  the  cross 
coupling  and  time  dependence  lost  by  linearization  can  be 
compensated  for  by  the  maneuver  covariance  matrix  Q.  Given 
the  target  acceleration  probability  model.  Figure  III-l,  the 
Q  matrix  can  be  calculated,  as  derived  in  Reference  3,  as 
time  varying  and  relates  the  cross  coupling  of  the  bearing 
and  range  channel  as: 

QR  =  f  0  0  0  and  QS  *  f  0  0  0  1  (55) 

000  000 

0  0  a2m  0  0  o  2m 

L  J  L  r  J 

where  o2«  -  acceleration  variance 

As  discussed  in  the  reference  the  Q(3,3)  element  can  be 
increased  to  make  the  missile  gain  matrix  put  more  weight  on 
any  target  acceleration  elements. 

If  the  ♦  matrix  is  maintained  as  time  varying  and 
nonlinear  then  the  Q  matrix  can  be  constant.  The  constant  Q 
matrix  can  be  calculated  as: 

Q  =  K  I  (56) 

j 

where  K  «  matrix  gain 

I  =  identity  matrix 

Since  the  reference  deals  extensively  with  the  time 
invariant  4>  and  the  time  varying  Q  matrix,  this  thesis  will 
deal  with  the  time  varying  *  and  constant  Q. 

If  Figure  III-4  is  used  to  define  the  maneuver 
probability  then  a  non-zero  mean  is  established.  The  Q 
matrix  maintains  the  same  properties  just  discussed  with  a 
different  calculation  for  o2«.  The  non-zero  mean  can  be 
implemented  by  increasing  Q(3,3),  increasing  the  weighting 
matrix  5,  or  by  not  assuming  the  U(k)  term  is  zero. 
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Kalman 

A 

xc  |k ) 

X(k*»  |k 

A 

Y(k«  »  |k 
PC  |k-i 
G  (k) 

PC  |k  ) 


filter  equations  for  this  third  order  model 

A  A 

«  XC  |k-i  )  +  G(k  )  •  [  Y(k)  -  YC  |k-i  )  ] 

A  A 

)  *  *.xc  |k )  +  r-u(k) 

)  =  H  •  X  C  *  1  |k  ) 

)  —  $*PC|k-i)*$T  +  5  *Q*  5T 

*  PC  |k-i  )  -HT  •  [  H-PC  |k-i  )  *HT  +  R(k)]  1 
=  P(k  |k-i  >  -  G  (k)  •  H* P  C  |  k  -  i  ) 


are : 

(57) 

(58) 

(59) 

(60) 
(61) 
(62) 
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VI.  SIMULATION 


To  aid  in  the  efforts  of  simulation,  the  Dynamic 
Simulation  Language  (DSL)  was  used  to  integrate  the 
equations  of  motion  for  the  target  and  missiles,  as  well  as 
the  antenna  angle  channel.  Two  programs  were  written  and 
are  listed  in  Appendix  B.  The  first  program  is  the  time 
varying  third  order  model.  The  second  program  is  the  time 
invariant  third  order  model  and  the  second  order  model. 

The  output  of  the  simulation  is  a  set  of  graphs  to 
compare  the  three  missile  models  and  their  effectiveness  in 
tracking  the  target. 

The  Kalman  Filter  is  implemented  in  a  Fortran  subroutine 
at  the  end  of  the  DSL  main  program.  The  basic  filter 
equations  used  were  described  in  the  previous  section. 

A.  ASSUMPTIONS 

The  following  assumptions  are  made  to  simplify  the 
implementation  of  the  Kalman  Filter  and  determine  the 
effects  of  the  time  varying  third  order  model. 

-  "Ideal"  missile  autopilot. 

-  Inertial  cartesian  reference  frame  for  angle 
measurements . 

-  Final  portion  of  intercept  only. 

-  Cross  coupled  effects  of  missile  motion  on  antenna 
stabilization  system  disregarded. 

-  Missile  initialized  to  collision  course. 

-  Missile  constant  speed  of  1500  kts  (Mach  1.5)  or 
2500  ft/sec. 

-  Target  constant  speed  of  400  kts  (Mach  .75)  or  667 
ft/sec . 

-  Target  lateral  acceleration  applied  perpendicular 
to  target  velocity  vector. 

-  Angle  of  Attack  not  accounted  for  in  velocity 
vector  calculations. 

-  Missile  located  at  the  center  of  the  cartesian 
reference  frame. 


-  Target  located  to  the  right  of  the  origin  of  the 
reference  frame. 

-  No  noise. 

B.  INITIALIZATION 

The  user  is  asked  at  the  beginning  of  the  program  to 
establish  the  geometry  by  giving  the  target  initial 
position,  heading,  speed  and  acceleration.  The  target 
heading  is  oriented  relative  to  a  vertical  line,  parallel  to 
the  Y  axis,  defining  the  North  or  000  heading.  The  initial 
missile  headings  are  calculated  for  constant  velocity,  zero 
acceleration  collision  course  with  a  time  to  go  estimate  of 
range/missile  velocity.  Antenna  parameters  are  initialized 
to  initial  line  of  sight  and  zero  angular  rate. 

C.  SECOND  ORDER  MODEL 

The  second  order  model  is  implemented  using  a 
proportional  navigation  constant  of  four  and  two  s-plane 
poles  at  s=-10.  This  gives  the  equation  for  angular 
acceleration  as: 

p'=  -20  3  +  100  (LOS  -  3)  (63) 

where  3  =  the  antenna  angular  acceleration 

3  =  the  antenna  angle  rate 

3  =  the  antenna  angle  position 

LOS  *  the  actual  angle  to  the  target 

D.  THIRD  ORDER  MODEL 

The  Kalman  filter  is  used  to  implement  the  third  order 
model.  The  DSL  main  program  calls  the  Kalman  Filter 
subroutine  at  a  sampling  time  of  h=0.01  seconds.  The  Kalman 
Filter  is  executed,  then  control  is  passed  back  to  the  DSL 
program.  A  proportional  navigation  constant  of  four  is  used 
as  discussed  in  Section  II. 
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1 .  Time  Invariant,  Constant  Phi  Model 

The  discrete  time  invariant  third  order  model 
divided  into  two  Kalman  filters  for  range  and  bearing  is: 


RNG  ( k  4  1  |  k  )  -  RPHI- RNG  (k  |k  ) 

(64) 

RNG  (k  |k 

)  *  RNG  (k  |k-i  >  +  GR(k  )  •  [  DELR  J 

(65) 

S(k4‘  |k 

)  =  SPHI  •  S  (k  |k  ) 

(66) 

S(“  |k  ) 

*  S  (k  |k-  l  )  +  GS(k  )  •  [  SDEL  ] 

(67) 

where 

RPHI  * 

SPHI  =  f  1  -01  .0005  ] 

0  1  .01 

0  0  1 


DELR  =  RM  -  RKP1 

RDM  -  RDKP1 

SDEL  =  f  LOS  -  SKP1  ] 

I  LOSD  -  SDKP1  I 

Using  Matlab  functions  of  Aker  and  Place,  with 
eigenvalues  of  0.5, 0.5  and  0.5  constant  gain  matrices  were 
obtained  for  range  and  bearing,  given  by  GR  and  GS, 
respectively . 


GR  = 


0.5  0.0125 

0.0025  1.0 
0.125  24.9 


GS 


1.5 

12.5 

1250.0 


2 .  Time  Variant.  Variable  Phi  Model 

The  discrete,  time  variant  third  order  model  divided 
into  two  Kalman  filters  for  range  and  bearing  is: 

RNG  (k  |k  )  =  RNG  (k  |  k  -  i  )  +  GR(k  )  •  [  DELR  J  (68) 

RNG  ( k  4  1  | k  )  =  RPHI  »RNG  (k  jk  ) 


(69) 


I 


PR  (k  |k  -  i  )  =  RPHI*PR(k  |  k  -  1  )  *RPHIT  +  QR 

GR (k)  *  PR ( k  | k  -  »  ) • HRT  £  j HR •  PR  ( k  j k -  1  ) • HRT  +RMCOV ( kl  | 1 

PR  (k  |k  )  =PR(k|k-i)  -  GR  (k)  *HR  •  PR  (k  jk  -  i  ) 

S(k|k)  «  S(k|k-i)  +  GS(k)*[  SDEL  ] 

S(ktl  |k  >  =  SPHI  •  S  {k  |k  ) 

PS(k|k-i)  =  SPHI*PS(k  |k-i  )  *SPHIT  +  QS 


(70) 


(71) 

(72) 


(73) 

(74) 

(75) 


GS  (k )  -  PS<“  |k-i  )  *HST  •  ^HS*PS(k  |k-i  )  ’•HST+RSCOV(k)J~1 
PS(“|k)  *  PS(k[k-i)  -  GS  (k)  -HS  •  PS  (k  |k-  i  ) 


(76) 

(77) 


The  time  variant  model  uses  two  different  4> 
matrices,  one  for  range  and  the  other  for  bearing.  Two 

other  matrices  must  be  specified  for  each  filter,  the 
initial  error  covariance  matrix  and  the  target  maneuvering 
covariance  matrix.  The  two  #  matrices  are: 


RPHI 


1  .01  .00005 

. 00005*A  1+ . 00005  *B  .01 

. 01  * A  . 01*B  1+ . 00005*B 


(78) 


SPHI  = 


where 


1  .01 

0  1+.00005*C 

0  . 01*C+ . 00005*D 


.00005 

. 01+ . 00005  *D 
1+.01*D+. 00005* (C+D) 


(79) 


A 

B 


■38  (8  ) 
■38  2 


C  =  -3r  /r 
D  =  -3r/  r 


The  initial  error  covariance  matrices  are  given  as: 


PR  (°  | o  ) 


500  0  0 

0  500  0 

0  0  500 


i 


i 


i 


i 
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PS  (°  |0  ) 


le4  0  0 

0  le4  0 
0  0  le4 


The  maneuvering  covariance  matrix  can  either  be  time 
varying  or  constant  as  discussed  previously.  The 
maneuvering  covariance  matrix  accounts  for  the  capabilities 
of  the  target  as  discussed  in  Section  IV.  Reference  3  gives 
the  derivation  for  the  time  varying  matrix  for  a  constant  4> 
matrix.  The  $  matrix  also  accounts  for  any  time  variance  of 
the  target  model  so  the  maneuvering  covariance  matrix  can  be 
constant.  A  constant  matrix  is  assumed  since  4>  is  time 
varying.  The  maneuvering  covariance  matrix  (QR  AND  QS)  are 
given  as: 


QR 


500  0  0 

0  500  0 

0  0  500 


QS 


.01  0  0 
0  .01  0 
0  0  .01 


The  simulation  calculates  the  range  model  then  the 
bearing  model.  The  results  of  each  filter  are  used  in  the 
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other  filter  to  calculate  the  $  values. 


E.  RESULTS 

Similar  simulations  of  the  ideal  missile  cases.  Section 
II,  were  used  to  evaluate  the  missiles.  The  simulation 
consists  of  head-on,  tail  and  beam  aspects  with  OG,  3G  and 
6G  target  acceleration.  The  results  of  the  computer  simu¬ 
lation  is  shown  in  graph  form  in  Figure  VI-1  through  VI-31. 
The  computer  program  listings  are  contained  in  Appendix  B. 

1 •  Gains 

Gain  comparison  plots  are  given  in  Figure  VI-1 
through  Figure  VI-4.  The  resultant  gains  from  the  time 
variant,  varying  Phi  third  order  missile  are  the  same  for 
each  scenario. 

In  predicting  R(k41  jk > ,  the  varying  Phi  model 
weights  the  range  error  by  .5  and  the  range  rate  error  by  0. 
The  constant  Phi  model  uses  weights  of  .5  and  .0125, 
approximately  the  same.  Figure  VI-1. 

Gains  for  predicting  RD{k41  |k  )  ,  for  varying  Phi,  are 
.001  and  .6  while  those  of  the  constant  Phi  are  .0025  and 
1.0,  Figure  VI-2.  Little  emphasis  is  placed  on  the  range 
error,  because  the  radar  is  measuring  range  rate,  with  a 
higher  weighting  factor. 

In  predicting  RDD(k  +  1  |k ) ,  small  gains  are  calculated 
by  the  varying  Phi  while  the  constant  Phi  model  places  a 
high  emphasis  on  range  rate  error.  The  noise  of  the  system 
will  be  noticed  more  in  the  prediction  of  RDD(k41  |k )  than 
the  other  parameters,  because  of  the  higher  weighting 
factor . 

Bearing  channel  gains  give  unusual  curves  for  the 
varying  Phi  model,  Figure  VI-4.  There  is  little  weight 
placed  on  non-observed  parameters  as  in  SG2  and  SG3 ,  during 
most  of  the  intercept,  except  in  the  initial  stage  and  final 
stage.  The  gains  are  highest  during  the  critical  stages  of 
the  intercept.  SGI  is  a  constant  1.0  giving  equal  weighting 
to  the  current  estimate  and  the  error.  The  constant  Phi 
model  has  higher  gains  giving  more  weight  to  any  errors. 
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Figure-  VI-4  Bearing  Covariance  Gains  SGI,  SG2 
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Figure  VI-5  Position  Flot  for  Head-on 
Aspect  No  Target  Turn 
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Figure  VI-8  Position  Plot  for  Head-on 
Aspect  3G  Target  Turn 
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REPRODUCED  AT  GOVERNMENT  EXPENSE 


Figure 
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REPRODUCED  AT  GOVERNMENT  EXPENSE 


Figure  VI-12 


Line  of  Siaht  Flot  for  Head-on 
Aspect  6G  Target  Turn 
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Fiaure  VI-13  Missile  Commanded  Acceleration  Plot  for  Head-on 

Aspect  3G  Target  Turn 
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Fiaure  VI-15  Line  of  Sight  Anale  Plot  for  Tail 
Aspect  No  Target  Turn 
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Figure  VI-23  Position  Plot  for  Beam 
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2 .  Head-on  Aspect 

Figures  VI-5  through  VI-16  are  the  results  of 
missile  guidance  comparisons  for  head-  on  aspect  initial 
condition  with  OG,  3G  and  6G  constant  target  acceleration. 
The  missile  begins  at  the  origin  of  the  graph.  The  target 
initial  position  is  at  x  =  10000  ft,  y  *  1000  ft.  Applied 
lateral  target  acceleration  is  away  from  the  missile. 

a.  0G  Target  Acceleration 

Flight  paths  for  the  three  models  are  shown  in 
Figure  VI-5.  All  three  paths  appear  to  be  the  same,  within 
the  accuracy  of  the  plotter.  The  line  of  sight  angle, 
Figure  VI-6  remains  relatively  constant  for  all  three 
models,  at  the  initial  value,  giving  a  constant  bearing 
decreasing  range,  no  manuever  intercept. 

Commanded  Missile  Acceleration,  Figure  VI-7, 
shows  that  the  second  order  model  pulls  more  lateral 
acceleration  than  third  order  models.  Higher  gain  terms 
make  the  constant  Phi  third  order  model  erratic  when 
compensating  for  initialization  error. 

b.  3G  Target  Acceleration 

In  the  position  plot,  Figure  VI-8,  the  second 
order  model  begins  to  lag  the  third  order  model.  The 
lagging  means  the  missile  is  slower  to  compensate  for  line 
of  sight  rates.  This  is  further  illustrated  by  Figure  VI-9, 
the  line  of  sight  angle  plot.  The  change  in  line  of  sight 
is  greater  for  the  second  order  model. 

Commanded  acceleration.  Figure  VI-10  shows  a 
larger  increase  for  the  second  order  model.  The  second 
order  model  requires  approximately  7 . 5G  to  intercept  a  3G 
target  while  the  third  order  model  require  4.66G. 
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c.  6G  Target  Acceleration 

Increased  target  acceleration  increases  the  lag 
of  the  second  order  model,  Figure  VI-11.  Magnitude  of  line 
of  sight  is  larger  for  the  second  order  model.  Figure  VI-12, 
but  the  missile  is  compensating  for  the  errors.  The  two 
third  order  models  are  fairly  close.  The  difference  is  only 
fractions  of  radians. 

The  constant  Phi  third  order  model  has  less 
initial  oscillations  in  commanded  acceleration  when  the 
target  acceleration  increases.  The  second  order  model 
requires  approximately  11G  for  a  6G  target,  7 . 5G  is  required 
for  the  third  order  models. 

3 •  Tail  Aspect 

Figures  VI-13  through  Figure  VI-21  are  the  results 
of  missile  guidance  comparisons  for  tail  aspect  initial 
condition  with  OG,  3G  and  6G  constant  target  acceler at ior. . 
The  missile  begins  at  the  origin  of  the  graph.  The  target 
initial  position  is  at  x  =  15000  ft,  y  =  -500  ft.  Applied 
lateral  target  acceleration  is  into  the  missile. 

a.  0G  Target  Acceleration 

As  shown  in  Figure  VI-13,  there  is  little 
difference  in  flight  paths  for  the  three  missiles  in  the  no 
target  acceleration  case.  Line  of  sight  angles  remain 
constant.  Figure  VI-14,  throughout  the  intercept. 

There  is  some  slight  difference  in  commanded 
accelerations.  Figure  VI-15,  with  the  second  order  model 
being  the  smaller  of  the  three  models. 

b.  3G  Target  Acceleration 

As  the  target  applies  acceleration,  the  lag  of 
the  second  order  missile  gives  a  shorter  flight  path,  than 
the  third  order  missiles.  Figure  VI-16.  Line  of  sight 
angles  are  smaller  for  the  second  order  model.  Figure  VI-17. 

Acceleration  required  for  the  second  order  model 
is  3.4G  and  5 . 3G  for  the  third  order  models,  Figure  VI-18. 
By  lagging  the  other  missiles,  the  second  order  model  allows 
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the  target  to  complete  part  of  the  intercept,  which  allows 
the  missile  to  pull  less  acceleration. 

c.  6G  Target  Acceleration 

With  higher  target  acceleration,  the  flight 
paths.  Figure  VI-19  have  approximately  the  same  differences 
as  the  3G  case.  The  second  order  model  maintains  a  better 
flight  path  throughout  the  intercept.  The  slope  of  the  line 
of  sight  curves.  Figure  VI-20,  are  higher  with  the  second 
order  model  maintaining  a  smaller  angle  difference. 

Commanded  acceleration  is  much  smaller  for  the 
second  order  model  than  the  third  order  models.  Figure  VI- 
21.  To  intercept  the  6G  target,  the  second  order  model 
requires  6.2G  and  the  third  order  models  require  8.1G. 

4 .  Beam  Aspect 

Figures  VI-22  through  VI-30  are  the  results  of 
missile  guidance  comparisons  for  head-  on  aspect  initial 
condition  with  OG,  3G  and  6G  constant  target  acceleration. 
The  missile  begins  at  the  origin  of  the  graph.  The  target 
initial  position  is  at  x  *  15000  ft,  y  =  0  ft.  Applied 
lateral  target  acceleration  is  away  from  the  missile. 

a.  OG  Target  Acceleration 

Figure  VI-22  through  VI-24,  show  the  three 
missiles  are  practically  the  same  for  a  no  target 
acceleration,  beam  .aspect  intercept.  All  three  missiles 
maintain  a  constant  bearing  decreasing  range,  small 
acceleration  intercept. 

b.  3G  Target  Acceleration 

Only  slight  differences  are  noticed  when  target 
acceleration  is  applied.  Flight  paths.  Figure  VI-25,  shows 
very  little  deviation.  The  line  of  sight  angle  difference 
of  .01  rad  between  the  models  is  approximately  .575  . 
Acceleration  is  higher  for  the  second  order  missile  to  allow 
it  to  fly  the  same  path  as  the  third  order  models.  Figure 
VI-27. 
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c.  6G  Target  Acceleration 

Flight  paths  have  a  pronounced  difference  with  a 
higher  target  acceleration,  Figure  VX-28.  Line  of  sight 
angles,  Figure  VI-29,  have  larger  magnitudes  and  higher 
slopes.  Commanded  acceleration  increases  dramatically,  now 
25G  is  required  of  the  second  order  model  and  9 . 3G  for  the 
third  order  models. 


VII .  CONCLUSIONS 


In  analyzing  second  and  third  order  missile  tracking  and 

guidance  subsystems,  the  following  conclusions  are  made: 

Proportional  navigation  guidance  is  the  optimum 
method  for  missiles,  given  current  design 
tradeoffs . 

Target  modeling  is  very  difficult  and  requires 
the  analysis  or  many  factors.  Acceleration 
probabilities  make  modeling  easier,  but  the 
proper  acceleration  model  must  be  chosen. 

Cross  coupling  between  coordinate  reference  axis 
components  does  exist  and  gives  errors  if  not 
accounted  for  in  the  system  model. 

Kalman  filters  are  the  best  predictors  for 
airborne  missiles,  if  one  is  required. 

Complete  time  varying  third  order  models  give 
better  results  than  approximated  linear,  time 
invariant  third  order  models. 

Only  small  differences  are  noticed  in  parameter 
values  between  second  and  third  order  models. 
Higher  accelerations  are  required  for  the  second 
order  model. 

Second  order  missiles  are  better  than  third  order 
missiles  in  tail  aspect,  constant  acceleration 
intercepts . 

Implementation  of  a  Kalman  filter  requires 
considerable  amounts  of  computer  resources,  with 
limited  time  to  complete  the  calculations. 

Some  parameter  terms  are  of  the  approximate  order 
as  system  noise  or  non  significant  calculations. 

Some  recommendations  for  future  study  and  consideration 


are : 


A  study  of  miss  distance  analysis  for  second  and 
third  order  models. 

Analysis  of  the  effects  of  the  Q  and  P  matrix 
initialization . 

Analyze  the  target  acceleration  probability  model 
to  find  optimum  values  to  assign  to  the 
probability  model. 

Determine  missile  cross  couple  effect  of  heading 
changes  and  autopilot  torques  on  the  sensor 
subsystem. 

Add  noise  to  the  system  to  determine  the  effects 
of  the  high  gain  terms  on  a  noisy  system. 
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APPENDIX  A.  IDEAL  MISSILE  PROGRAM  LISTING 


C  MISSILE  PROGRAM  FOR  FLIGHT  PATH  COMPARISON  IN  THE 

C  THESIS,  THIS  PROGRAM  IS  DESIGNED  TO  FOLLOW  THREE  IDEAL 

C  MISSILES  USING  DIFFERENT  GUIDANCE  TECHNIQUES  TO  TARGET 

C  INTERCEPT.  THE  GUIDANCE  TECHNIQUES  ARE:  PROPORTIONAL 

C  NAVIGATION,  PURE  PURSUIT  AND  LEAD  PURSUIT.  PARAMETERS 

C  ARE  ALSO  OBTAINED  FOR  A  FOURTH  MISSILE  USING  A  DIRECT 

C  FLIGHT  PATH. 

INITIAL 

CONST 

G=32 . 2 , D2R= . 0175, K2F=1. 66667 , PITCH=2 .7,PI=3.14159, RM1=20000 
K=0 
NN=0 

MISSX0=0.0 
MISSY0=0.0 
VM  =  2500 . 0 

READ  (2,5)  VT , AT , THDG , TGTXO , TGTYO 
5  FORMAT  (F6.1,2X,F5.1,2X,F6.1,2(2X,F10.2) ) 

TGTV=VT*  K2F 
TGTA=-AT*G 
THDG=THDG*D2R 
TGTVXO=TGTV*SIN (THDG) 

TGTVYG=TGTV*CGS (THDG) 

C 

C  INITIALIZATION  OF  DIRECT  INTERCEPT  MISSILE 

DXT=MISSXO 
DYT=MISSYO 
C 

C  INITIALIZATION  OF  PROP  NAV  MISSILE  CONSTANT  VELOCITY, 

ZERO  ACCEL 

R0=( (TGTXO  -  MISSXO ) *  *2  +  (TGTYO  -  MISSY0)**2  )**.5 
TTG0=  RO/VK 

LOS  =  ATAN2 (TGTYO-MISSYO.TGTXO-MISSXO) 

PHDG=ATAN2 ( TGTYO+TGTVYO  *TTG0-MISSY0 , - 
TGTXO+TGTVXO  *TTG0 -MISSXO ) 

BDO  =  TGTV*COS (THDG) / (COS (LOS)  *  RO) 

BO  =  LOS 
C 

C  INITIALIZATION  OF  LEAD  PURSUIT  MISSILE 

GO  =  BG 
GDO  =  0 
LPHDGO  =  LOS 

C  ****************** 

METHOD  RKSFX 

DERIVATIVE 

C 

C  TARGET  POSITION  UPDATING 

C  TGTHDG=INTGRL (THDG, (-1*AT) *PITCH*D2R) 

TGTAX=TGTA*COS (TGTHDG) 

TGTAY=-TGTA*  SIN (TGTHDG ) 

TVELX=INTGRL (TGTVXO , TGTAX) 

TVELY= INTGRL ( TGTVYO , TGTAV ) 

XT=INTGRL (TGTXO ,TVELX) 

YT=INTGRL (TGTYO ,TVELY) 

TGTHDG  =  AT AN 2 (TVELX , TVELY) 

C 

C  PROP  NAV  MISSILE  POSITION  UPDATING 

BDDOT  =  -20*BDOT  +  100  * ( PNLOS-B ) 

BDOT  =  INTGRL (BDO , BDDOT) 

B  =  INTGRL (BO , BDOT) 

PNHDG  =  INTGRL (PHDG , 4 *BDOT) 

PXM=INTGRL (MISSXO , VM*COS (PNHDG) ) 

PYM=INTGRL (MI SSYO , VM* SIN (PNHDG ) ) 

C 

C  PURE  PURSUIT  MISSILE 
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PPHDG  =  ATAN2 (YT-PPYM , XT-PPXM) 

PFXM  =  INTGRL  MISSXO,VM*COS (PPHDG) ) 

PPYM  =  INTGRL (MISSYO , VM* SIN (PPHDG) ) 

C 

C  LEAD  PURSUIT  MISSILE 

LPHDG 

AT AN 2 (YT+TVELY*0 . 5*TTG-LPYM , XT+TVELX*  0 . 5*TTG-LPXM) 

LPXM  =  INTGRL (MIS SXO , VM*COS (LPHDG) ) 

LPYM  =  INTGRL {MISSYO, VM* SIN (LPHDG) ) 

C 

Q  ******************* 

C 

DYNAMIC 

C 

C  PROP  NAV  MISSILE  GEOMETRY  UPDATE 

PNAM  =  4*BDOT*PNRD 

PNR=( (XT-PXM) **2  +  ( YT-PYM) *  *  2 ) *  * . 5 
PNLOS  =  ATAN2 (YT-PYM. XT-PXM) 

PNRD=TVELX*COS (PNLOS)  +VM*COS ( PNHDG-PNLOS ) 

PNLOSD  =  { -TVELX  *SIN ( PNLOS ) -VM*  SIN ( PNHDG-PNLOS ) ) /PNR 
C 

C  PURE  PURSUIT  GEOMETRY  UPDATE 

PPR  =  ( (XT-PPXM) **2  +  (YT-PPYM) *  *  2 ) *  * . 5 
PPLOS  =  ATAN2 (YT-PPYM, XT-PPXM) 

PPRD  =  TVELX*COS (PPLOS)  -  VM*COS ( PPHDG-PPLOS ) 

PPLOSD  =  -TVELX*SIN (PPLOS) /PPR 
C 

C  LEAD  PURSUIT  GEOMETRY  UPDATE 

LPR  =  (  (XT-LPXM) **2  +  (YT-  LPYM) *  *  2 ) *  *  .  5 
LPLOS  =  ATAN2 (YT-LPYM, XT-LPXM) 

LPRD  =  TVELX*COS (LPLOS)  -  VM*COS ( LPHDG-LPLOS ) 

LPLOSD  =  ( -TVELX*  SIN ( LPLOS ) 

VM*SIN (LPHDG-LPLOS) ) /LPR 

C 

c 


c 


50 

51 

52 
54 

C 

C 

C 

C 


C 
C 

SAMPLE 

C  SAVE  (A)  0.1 , XT, YT,PXM,PYM,PPXM, PPYM, LPXM, LPYM 
C  PRINT  1. 0, XT, YT, LPXM, LPYM, GDOT, LPLOS, LPHDG 
CONTROL  FINTIM=  8.0,DELT=.01 
TERMINAL 

DHDG  =  AT AN 2 (DYT-MISSYO , DXT-MISSXO ) 

C  GRAPH  ( A/A , DE=TEK61 8 )  XT 

( SC  =  1 600 , LO=0 .00)  , YT ( SC  =  7  50 , LO=0 .0,PO=16000) 

C  GRAPH  ( A / A , O V ) 

PXM (SC=1600 , LO=0 . 0 , AX=OMIT) , PYM ( SC=750 , LO=0 ) 


TTG  =-  LPR/LFRD 

DIRECT  INTERCEPT  MISSILE  GEOMETRY  UPDATE 
RD= (  (YT -MISSYO ) *  *  2  +  (XT-MISSX0 i *  *  2 ) *  * . 5 
DXT , DYT , FON=  CHECK ( RD , TIME , VM , XT , YT ) 

IF  (PXM  .GT.  XT)  CALL  ENDRUN 

IF  (K  .LE.  0.0)  THEN 
NN=NN+1 

WRITE  (31,50)  XT, YT, PXM, PYM 
WRITE  32,50)  PPXM , PPYM , LPXM , LPYM 
WRITE  (33,51)  TIME, PNLOS, PPLOS, LPLOS 
WRITE  34,52)  TIME , PNLOSD , PPLOSD , LPLOSD 
WRITE  (36,54)  TIME , PNHDG , PPHDG , LPHDG 
FORMAT (4 (F10.2.2X) ) 

FORMAT  (F5 . 2 , 3 ( 2X , F10 . 5  5  ) 

FORMAT  (F5.2. 3(2X,F10.5) ) 

FORMAT (F5.3,2X,3(F10.6,2X) ) 

K=10 

IF  ( PNR . LT • 0 . 1 *R0 )  THEN 
K=1 
ENDIF 
END  IF 
K=K-1 
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C  GRAPH  ( A/A , OV ) 

PPXM { SC=1600 . LO=G . 0 , AX=OMIT)  , PPYM ( SC  =  7  50 , LO=0 , AX=OMIT ) 

C  GRAPH  (A/A.OV) 

LPXM ( SC=1600 , LO=0 . 0 , AX=OMIT) , LPYM ( SC=750 , LO=0 , AX=OMIT) 

WRITE  (2,15)  DHDG , FON 

15  FORMAT  (F10.7 ,2X,F5.2) 

WRITE  (1,16)  NN 

16  FORMAT  (F4.0) 

END 

STOP 

FORTRAN 

SUBROUTINE  CHECK  (RD , TIME , VM , XT , YT , DXT , DYT , FON ) 

REAL* 8  RD , TIME , VM , XT , YT , DXT , DYT , FON , DMISS 
DMISS  =  VM*TIME 
IF  (DMISS  .LT.  RD)  THEN 
DXT  =  XT 
DYT  =  YT 
FON  =  TIME 
ENDIF 
RETURN 
END 
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C  THIS  IS  A  DSL  PROGRAM  TO  FIND  THE  LOS  AND  RANGE  OF  A  « 

C  TARGET  IN  A  CONSTANT  G  TURN  FROM  A  MISSILE  ON  A 

C  DIRECTPATH . 

C 

c 

INITIAL 

CONST 

G-32.2,D2R-.0175.PITCH-2.7 ,K2F«1. 66667 ,1-1.0, PI=3. 14159, K-0 

READ  <2,10)  VT , AT , THDG , TGTXO , TGTY  0  . 

READ  (2,11)  MISHDG , DONE 

10  FORMAT  (F6.1,2X,F5.2,2X,F6.2,2(2X.F10.2)) 

11  FORMAT  ( F10 . 7 , 2X , F5 . 2 ) 

C  MISSILE  PARAMETERS 

MISSXO  =0.0 
MISSYO  =0.0 
VM  =  2500.0 

C  TARGET  PARAMETERS  , 

TGTHDG=THDG*D2R 
TGTV=VT*K2F 
TGTG  =-AT*G 

TGTTVELXO-TGTV*  SIN (TGTHDG ) 

TGTTVELYO=TGTV  *  COS (TGTHDG) 

C  EQUATIONS  OF  MOTION 

METHOD  RKSFX 

DERIVATIVE 

TGTHDG=INTGRL (THDG, -AT*PITCH*D2R)  | 

TGTAX=TGTG*COS (TGTHDG) 

TGTAY=-TGTG*SIN (TGTHDG) 

TVELX=INTGRL (TGTTVELXQ , TGTAX ) 

TVELY=INTGRL ( TGTTVELYO , TGTAY ) 

XT=INTGRL (TGTXO , TVELX) 

YT=INTGRL (TGTYO , TVELY) 

C  MISSILE  POSITION  UPDATE 

XK  =  INTGRL (MISSXO, VM*COS (MISHDG) )  I 

YM  -  INTGRL (MI SSYO, VM* SIN (MI SHDG ) ) 

DYNAMIC 

R= •  • XT-XM) *  *  2  +  (YT-YM)  **2) **.5 
LOS  =  ATAN2 (YT-YM, XT-XM) 

RD  =  TVELX*COS (LOS)  -  VM*COS (MISHDG-LOS ) 

L OS!  =  (  TVELY*COS (LOS)  -  VM* SIN (MISHDG-LOS )) /R 
IF  iK.EQ.10)  THEN 
MM  =  MM+1 

WRITE  (37,15)  XM , YM 
WRITE  (38,16)  TIME, LOS ,LOSD, MISHDG 
15  FOR:. AT  (2(F9.2,3X;  ) 

It  FOR  .AT  (F5.2,2X,2(F10.5,2X)  ,F7.4) 

K=C 
END  IF 
K=R- 1 

IF  :TIME  .GT.  DONE)  CALL  ENDRUN 

SAMPLE 

CONTROL  FINTIM=10 . 0 , DELT= . 01 
C  PRINT  1.0, XT , YT , XM , YM , LOS , R 
C  SAVE  (D)  0.1, XT , YT , XM , YM , R , LOS 
TERMINAL 

REAI  (1.26)  NN 
WRITE  (1,26)  MM 
26  FOR'*  AT  (F5.0) 

C  GRAPH  '  (A/D, DE=TEK618)  XT ( SC  =  1 500 , LO=0 )  , YT ( PO=l 5000 ) 

C  GRAPH  ( A/D , DE=TEK6 1 8 , OV )  XM ( SC=1500 , LO=0 , AX=OMIT) , YM 
C  GRAPH  (A/D,DE=TEK618)  TIME, LOS 
C  GRAPH  ( A/D , DE=TEK618 )  TIME , R 
END 
STOP 
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APPENDIX  B.  THIRD  ORDER  SIMULATION  PROGRAM  LISTING 


2) , RNG ( 3 ) ,S(3) , DELR ( 3 ) 


DIMENSION 
SG  ( 3 ) 


C  PROP  NAV  MISSILE  PROGRAM  FOR  THESIS 

INITIAL 
D 

PS  ( 3 , 3 )  ,  PR  ( 3 , 3 )  , RMCOV ( 2 
D  DIMENSION  RG(3,3) 

K=0 

MM=0 

HH=0 

C 

METHOD  RKSFX 

CONST  G=32 .  2,D2R= . 0175 , K2F=1 .66667 
TK  =  .01 
MISSX0=0 . 0 
MISSY0=0.0 
VM  =  2500.0 
AMO  =0.0 

READ  (2,10)  VT , AT , THDG , TGTXO , TGTYO 
10  FORMAT  (F6.1,2X,F5.1,2X,F6.1,2X,F10.2,2X,F10.2! 

TGTV=VT*  K2F 
TGTA= ( -1  * AT ) *G 
THDG=THDG*D2R 
TGTVXO=TGTV* SIN (THDG) 

TGTVYO=TGTV*COS (THDG) 


C 

c 


INITIAL 
PS(1, 
PS(1, 
PSU, 
PS  (2 
PS  (2 
PS  (2 
PS  (3 


PS (0/0)  MATRIX 

1 )  =1 . OE+4 

2 )  =0 . 0 
3  )  =0 . 0 

1) =PS (1,2) 

2 )  =1 . OE+4 
3  )  =0 . 0 
1 ) =PS (1,3! 


C 

C 


c 

c 


INITIAL  PR (0/0)  MATRIX 
PR (1,1)  =  500 
PR (1,2)  =  0 
PR (1,3)  =  0 
PR (2,1)  =  PR ( 1 , 

PR ( 2 , 2  >  =  500 
PR (2, 3)  =  0.0 
PR (3,1)  =  PR ( 1 , 

PR (3, 2)  =  PR (2, 

PR ( 3 , 3 )  =  500 


C 

C 

MATRIX 


INITILIZE  THE 
RMCOV (1,1) 
RMCOV  1,2 
RMCOV (2,1) 
RMCOV (2, 2) 


2) 


3 

3) 


RANGE  MEASUREMENT  COVARIANCE  MATRIX 
=  0.0 
=  0.0 
=  RMCOV (1 
=  0.0 


2) 


C 

c 

ZERO 


INITIALIZE  THE  BEARING  MEASUREMENT  NOISE  COVARIANCE 
SMCOV  =0.0 


INITIALIZATION  OF  PROP  NAV  MISSILE  CONSTANT  VELOCITY, 
ACCEL 

LOS  =  ATAN2 (TGTYO  -  MISSYO  ,  TGTXO  -  MISSXO) 

R=  (  ( TGTXO  -  MISSXO )  *  *  2  +  (TGTYO  -  MISSY0)**2  )“.5 
TTG0=  R/VM 

PHDG=ATAN2 ( TGTYO +TGTVY0*TTG0 -MISSYO , TGTX0+- 
TGTVX0*TTG0 -MISSXO ) 

VMXO  =  VM*COS ( PHDG ) 
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VMYO  =  VM*SIN (PHDG) 

RKFI  =  R 

RDKF1  =  -VM*COS {PHDG-LOS )  +  TGTV* SIN (THDG ) /COS (LOS ) 

RDDKP1  =  0 
SKP1  =  LOS 

SDKP1  =  (TGTVYO/COS (LOS)  -  VM* SIN (PHDG-LOS )) /R 
SDDKP1  =  0 
BO  =  LOS 
BDO  =  0 
C 

RNG(l)  =  RKP1 
RNG ( 2 )  =  RDKP1 
RNG (3)  =  RDDKP1 
S(l)  =  SKP1 
S ( 2 )  =  SDKP1 

S ( 3 )  =  SDDKP1 

C 

DERIVATIVE 

C 

C  TARGET  POSITION  UPDATING 

TGTHDG  =  ATAN2 ( TVELX , TVELY ) 

TGTAX  =  TGTA* COS (TGTHDG) 

TGTAY= ( -1 *TGTA) *  SIN (TGTHDG) 

TVELX  *  INTGRL(TGTVXO, TGTAX) 

TVELY=INTGRL (TGTVYO , TGTAY ) 

XT  =  INTGRL (TGTXO , TVELX) 

YT=INTGRL (TGTYO , TVELY ) 

C 

C 

C  THIRD  ORDER  PROP  NAV  MISSILE  POSITION  UPDATING 

BDOT  =  INTGRL (BDO, BDDOT) 

6  =  INTGRL (BO, BDOT) 

AM  =  U 

MVELX  =  INTGRL (VMXO, -AM* SIN ( PNHDG ) ) 

MVELY  *  INTGRL (VMYO, AM*COS (PNHDG) ) 

PNHDG  =  AT AN 2 (MVELY, MVELX) 

PXM= INTGRL (MISSXC , MVELX ) 

F YM= INTGRL ( MI SSYO, MVELY) 

C 

DYNAMIC 

C 

C  THIRD  ORDER  PROP  NAV  MISSILE  GEOMETRY  UPDATE 

RM  =  ( (XT-PXM) *  *  2  +  ( YT-PYM) *  *  2 ) *  * . 5 
LOS  =  ATAK2 (YT-PYM, XT-PXM) 

RDOTM  =  TGTV*SIN (TGTHDG) /COS (LOS) 

VM*COS (PNHDG-LCS) 

RDDOTM  =  TGTA*COS (TGTHDG-LOS) 

LOSD  =  (  TVELX* SIN ( LOS )  +  VM* SIN ( PNHDG-LOS ) ) /RM 
LOSDD  =  TGTAY*COS (LOS) /RM 
C  COMPUTE  THE  ERROR  TERMS 

DELR(l)  =  RM  -  RKP1 
DELR ( 2 )  =  RDOTM  -  RDKP1 
SDEL  =  LOS  -  SKP1 
C 

c 

CALL 

KALMAN ( RNG, RM , RDOTM, RDDOTM, RK , RDK , RDDK , RKP1 , RDKP1 , RDDKP1 , . . . 


R, DELR, PR, RMCOV,S, LOS, LOSD, LOSDD, SK, SDK, SDDK, TIME, SKP1, . . . 
SDKP1 , SDDKP1 , SDEL , SMCOV , P  S , TK , RG , SG , HH ) 
BDDOT=SDDKP1+10* ( SDK-BDOT) +33 . 33333* (LOS-B) 

U  =  -  4  *  BDOT  *  RDOTM 
C 

C  SAVE  VALUES  OF  THE  GAIN  MATRICES 

GR11=RG (1,1) 

GR1 2=RG ( 1 , 2 
GR21=RG (2,1) 

GR22=RG (2,2) 
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GR31=RG  (3,1) 

GR32=RG (3,2) 

GS1=SG ( 1 ) 

GS2=SG(2) 

GS3*SG(3) 

IF  (  PXM  .GT.  XT  )  THEN 
CALL  ENDRUN 
END  IF 
C 

IF  (K  .LE.  0)  THEN 
MM=MM+1 

WRITE  (40,20)  XT , YT , PXM , PYM 

WRITE  42,21  TIME ,  GR11 ,  GR12  ,  GR21 ,  GF.22  ,  GR31 ,  GR32 
WRITE  (44,22)  TIME , GS1 , GS2 , GS3 
WRITE  (46,23)  TIME , LOS , LOSD , LOSDD , U 
K=10 

IF  (RM  .LT.  .  1*R)  THEN 
K=3 
END  IF 
END  IF 
K=K- 1 

20  FORMAT (4 (2X.F10.2) ) 

21  FORMAT  (F5.2,6(1X,F10.4) } 

22  FORMAT  (F5 . 2 , 3 (2X, E12 . 5) } 

23  FORMAT  (F5 . 2 , 3 ( 2X , Ell . 4 ) , 2X , E14 . 6 ) 

SAMPLE 

C  STATEMENTS  TO  SAVE  DATA  FOR  USE  WITH  GRAFAEL 

C  SAVE  (A)  0.1, XT , YT , PXM , PYM 
C  SAVE  (B)  0.1,  LOS , B , SK 
C  SAVE  (C)  0.1, LOSD, BDOT, SDK 

C  SAVE  D  0.1, SDEL 

C  SAVE  (E;  0.1, RM , RK , RDDOTM , RDDK 
C  SAVE  (F)  0 . 1 , U , BDDOT 
C  SAVE  (G)  0.1, GS1 , GS2 , GS3 
C  SAVE  (H)  0.1, GR11 , GR12 , GR21 , GR22 

C  PRINT  0.1, RM , RKP1 , RK , RDOTM , RDKP1 , RDR , RDDOTM , RDDKP1 , RDDK 

C  PRINT  0.1.PNHDG, LOS, B,SK, LOSD, BDOT, SDK, LOSDD, BDDOT, SDDK 

CONTROL  FINTIM=10 . 0 , DELT= . 01 
TERMINAL 

WRITE  (1,30)  MM 
30  FORMAT  (F6.1) 

C  STATEMENTS  FOR  PLOTTING  WITH  GRAFAEL 

C  GRAPH  ( A/A , DE=TEK61S )  XT 

(SC=1600,LO=0 .0) ,YT(SC=500,PO=16000) 

C  GRAPH  ( A / A , O V )  PXM  ( SC=1600 , AX=OMIT) , PYM ( SC=500 ) 

C  GRAPH  (B/B,DE=TEK618)  TIME , LOS ( SC= . 025 , LO=- . 1 ) 

C  GRAPH  B/B,OV)  TIME (AX=OMIT) , B ( PO=7 . 5, SC= . 025 , LO=- . 1 ) 

C  GRAPH  (B/B.OV)  TIME ( AX=OMIT ) , SK ( AX=OMIT , SC= . 025 , LO=- . 1 ) 

C  GRAPH  (C/C , DE=TEK618 )  TIME, LOSD 
C  GRAPH  (C/C.OV)  TIME ( AX=OMIT ) , BDOT ( PO=7 . 5 ) 

C  GRAPH  C/C,OV)  TIME (AX-OMIT), SDK (AX=OMIT) 

C  GRAPH  (D/D , DE=TEK6 18 )  TIME, SDEL 
C  GRAPH  ( E/E , DE=TEK6 1 8 )  TIME , RM ( SC=2000 . 0 . LO=0 . 0 ) 

C  GRAPH  ( E / E , OV )  TIME  ( AX»OMIT )  , RK  ( SC  =  2OO0 . 0 , LO=0 . 0 ) 

C  GRAPH  (F/E ,  DE=TEK618  )  TIME  ,  F.DDOTM  (  SC*500 . 0 ) 

C  GRAPH  ( F / E , OV )  TIME (AX=OMIT) , RDDK (AX=OMIT) 

C  GRAPH  G/F,DE=TEK618)  TIME,U 

C  GRAPH  (H/G , DE=TEK61 8 )  TIME,GS1 

C  GRAPH  ( I/G , DE=TEK618 )  TIME , GS2 

C  GRAPH  ( J/G , DE=TEK618 )  TIME,GS3 

C  GRAPH  ( K/H , DE=TEK61 8 )  TIME , GR11 

C  GRAPH  ( L/H , DE=TEK6 1 8 )  TIME , GR12 

C  GRAPH  M/H,DE=TEK618  TIME , GR21 

C  GRAPH  (N/H , DE=TEK618 )  TIME , GR2  2 

C  GRAPH  (O/F , DE=TEK618 )  TIME, BDDOT 
END 
STOP 
FORTRAN 
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SUBROUTINE  KALMAN ( RNG , RM , RDOTM , RDDOTM , RK , RDK , RDDK , 

RKP1 , RDKP1 , RDDKP1 , K , DELR . PR , RMCOV , S , LOS , LOSD , LOSDD , 

SK , SDK , SDDK , TIME , SKP1 , SDKPl , SDDKP1 ,  SDEL , SMCOV , 

PS , TK , RG , SG , HH ) 

C  SUBROUTINE  TO  ITERATE  A  KALMAN  FILTER  FOR  RANGE 

C  VARIABLES 

C  GIVEN  THE  COVARIANCE  MATRIX  AND  OBSERVATIONS 

REAL *8  RNG 13) , RM.RDOTM , RDDOTM . RK , RDK , RDDK , RKP1 , 

*  RDKP1 , RDDKPl , DELR ( 3 )  , PR  (  3 , 3 )  , RMCOV ( 3 , 3 )  , S  (  3 )  , 

*  LOS , LOSD, LOSDD, SK, SDK, SDDK, SKP1 , SDKPl , SDDKP1 , SDEL , 

*  SMCOV, TIME,  PS  ( 3 , 3)  , TK, TKSQ , A , B , C , D , 

*  TEMPI (3), 
TEMP2 (2,2) , TEMP 3 (3,3) ,RPHI(3,3) ,COVR(3,3) , 

*  DET, SCOV.SPHI (3, 3) , RG (3 , 3) , SG (3) , QR ( 3 , 3) , QS ( 3 , 3) 

*  TKSQ  =  TK‘TK 
C 

C  FIND  THE  NEW  VALUES  OF  RPHI  FROM  THE  PREVIOUS  VALUES 

OF  SIGMA 
C  MATRIX 

A  =  -3*SDKP1*SDDKP1 
B  =  -3*SDKP1*SDKF1 
RPHI (1,1)  =  1 
RPHI (1 , 2)  =  TK 
RPHI  1,3)  =  . 5*TKSQ 
RPHI (2,1)  =  . 5*A*TKSQ 
RPHI (2,2)  =  1  +  .5*B*TRSQ 
RPHI (2, 3)  =  TK 
RPHI (3,1)  =  A*TK 
RPHI (3,2)  =  B*TK  +  .5*A*TKSQ 
RPHI (3,3)  -  1  +  . 5*B*TRSQ 
C 

C  FIND  THE  PROJECTED  COVARIANCE  PR(K/K-1)  = 

RPHI*PR(K-1/K-1) *RPHI 

CALL  ZERO {TEMP 3, 3) 


DO  104  L=1 , 3 
DO  103  M=1 , 3 
DO  102  N=1 . 3 

TEMP  3 ( L , M ; =TEMP  3 ( L , M )  +  RPHI (L,N) *PR (N ,M) 
CONTINUE 
CONTINUE 
CONTINUE 

CLEAR  THE  OLD  COVARIANCE  MATRIX 
CALL  ZERO (PR, 3) 

MULTIPLY  BY  RPHI  TRANSPOSE 
DO  107  L=1 , 3 
DO  106  M=1 , 3 
DO  105  N-l , 3 

PR ( L , M) =PR ( L , M ) +TEMP3 ( L , N) ‘RPHI  (M,N) 
CONTINUE 
CONTINUE 
CONTINUE 

DEFINE  THE  Q  MATRIX  OF  MANEUVER  COVARIANCE 


120 


C  NOW  PR  IS  THE  COVARIANCE  OF  THE  PREDICTED  ESTIMATE 

PRIK/K-1) 

C  FIND  THE  ESTIMATE  OF  RANGE  MATRIX  AT  STEP  K 

C 

C  ZERO  A  TEMPORARY  MATRIX 

CALL  ZERO  (TEMP 2, 2) 

DO  121  L=1 , 2 
DO  120  M=1 , 2 

TEMP 2 (L ,M)  =  PR (L ,  M)  +  RMCOV(L.M) 

120  CONTINUE 

121  CONTINUE 

DET  *  TEMP2 (1,1) *TEMP2 (2,2)  -  TEMP2 (1,2) *TEMP2 (2,1) 


COVR (1,1)  =  TEMP 2 (2,2) /DET 
COVR (1,2)  =  (-1) *TEMP2(1,2> 
COVR (2,1)  «=  ( -1 )  *TEMP2{2,1) 
COVR (2, 2)  »  TEMP2U.1) /DET 
HERE  COVR  =  (HPH  +  R)  INVERSE 
ZERO  A  TEMPORARY  MATRIX 
CALL  ZERO  (RG , 3) 

DO  132  L=1 , 3 
DO  131  M=1 , 2 
DO  130  N=1 , 2 

RG (L.K)-RG (L.M)  + 
CONTINUE 
CONTINUE 

RG (L , 3 )  =  0 
CONTINUE 

RG  =  PH (HPH  +  R)  ( 

ZERO  A  TEMP  MATRIX  FOR 
DO  140  L=1 , 3 
TEMPI ( L ) — 6 
CONTINUE 


/DET 

/DET 


PR (L , N )  *  COVR (N , M) 


R)  (INVERSED)  H 
FOR  THE  RNG  MATRIX 


DO  142  L-l , 3 
DO  141  M=1 , 2 

TEMPI (L)  =  TEMPI (L)  +  RG (L , M)  *  DELR(M) 

CONTINUE 
CONTINUE 
DO  143  N=1 , 3 

RNG(N)  =  TEKFl(N)  +  RNG (N) 

CONTINUE 

SAVE  THE  VALUES  OF  RANGE  MATRIX  AT  STEP  K 
RK  =  RNG ( 1 ) 

RDK  =  RNG ( 2 ) 

RDDK  =  RNG ( 3 ) 

ZERO  THE  OLD  RANGE  TEMPORARY  MATRIX 
DO  150  N  =  1.3 
TEMPI (N)  ■  6 
CONTINUE 

FIND  THE  ESTIMATE  OF  THE  STEP  K+l  FOR  THE  RANGE  MATRIX 


DO  152  L  = 
DO  151  M  = 
TEMPI (L)  * 
CONTINUE 
CONTINUE 


1.3 

1.3 

<  TEMPI (L) 


RPHI (L.M) 


RNG (M) 


SAVE  THE  VALUES  OF  RNG ( K+l/K) 

DO  153  N=1 , 3 

RNG(N)  =  TEMPI (N) 

CONTINUE 

RKP1  =  RNG ( 1 ) 

RDKP1  =  RNG (2) 

RDDKP1  «  RNG (3) 

FIND  THE  COVARIANCE  OF  FILTERED  ESTIMATE 
RG  =  PH (HPH+R) INVERSED  H 
THEREFORE  P(K/K)  =  P(K/K-1)  +  RG*P(K/K-1) 
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CALL  ZERO (TEMP 3, 3) 

DO  163  L=1 , 3 
DO  162  M=l,3 
DO  161  N=1 , 3 

TEMP 3 (L ,  M)  =  TEMP3 (L,M)  +  RG (L , N ) ‘PR (N , M) 
CONTINUE 
CONTINUE 
CONTINUE 

DO  165  L=1 , 3 
DO  164  M=l,3 

PR ( L , M)  =  PR (L ,M)  -  TEMP 3 (L,M) 

CONTINUE 

CONTINUE 

NOW  PR  IS  THE  COVARIANCE  OF  FILTERED  ESTIMATE  P(K/K) 


VARIABLES 


SUBROUTINE  TO  ITERATE 


KALMAN  FILTER  FOR  SIGMA 


GIVEN  THE  COVARIANCE  MATRIX,  OBSERVATION. 
WHERE  H=  (1  0  0) 

CALCULATE  THE  NEW  SPHI  MATRIX 
C  =  -3*RDDK/RK 
D  =  -3  *RDK/RK 
SPHI (1,1)  =  1 
SPHI (1,2)  =  TK 
SPHI (1,3)  =  . 5*TKSQ 
SPHI (2,1)  =  0 
SPHI (2,2)  *  1  +  . 5*C‘TKSQ 
SPHI  2,3  =  TK  +  . 5*D*TKSQ 

SPHI (3,1)  =  0 

SPHI  3,2  ■  C*TK  +  . 5‘C*D‘TKSQ 

SPHI (3,3)  «  1  +  D*TK  +  .5* (C+D) *TKSQ 


SPHI 


1  +  D*TK  +  . 5* (C+D) *TKSQ 


CALCULATE  THE  NEXT  PROJECTED  P  MATRIX  FOR  NEXT  STEP 

T 

PS  (K, K-1)=SPHI*PS*SPHI 
CLEA:  A  TEMPORARY  MATRIX 


2C  . 

P.  202  L=1 , 3 
DO  201  M=1 , 3 
DC  200  N=1 , 3 
TEMP 3 ( L , M) 
CONTINUE 

2i  _ 

C  C  NTINUE 

2C  . 

CTNTINUE 

C 

ZERO  THE  OLD  PS 

20  6 

C-r.LL  ZERO  (PS 
DO  205  L=1 , 3 
DO  204  M®=1 , 3 
DO  203  N-1.3 
PS (L,M) =PS 
CONTINUE 

204 

CONTINUE 

20  5 

CONTINUE 

DEFINE  THE  Q  MATRIX  OF  MANEUVER  COVARINCE 
OS  ( 1 , 1 )  =  0.01 
06  (1,2;  =  0.0 


OS  1,3)  *  0.0 
06(2,1)  =  OS (1,2) 
03(2,2)  =  0.01 
03(2,3)  =  0.0 

06(3,3)  =  0.01 
DO  211  L=1 , 3 
DO  210  M=1 , 3 


PS (L , Mj  =  PS (L,M)  +  QS (L,M) 
CONTINUE 
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C 

c 

c 


CONTINUE 

T 

CALCULATE  (HPH  +  R)  INVERSE 


c 

r 

SCOV  = 

PS (1,1)  +  SMCOV 

c 

SG  ( 1 ) 

=  PS(l,l)/SCOV 

SG  (2) 

=  PS(2,1  /SCOV 

SG  3 

=  PS (3,1  /SCOV 

c 

NOW  FIND  THE  CURRENT  VALUES  < 

S  (1)  = 

S ( 1 )  +  SG ( 1 ) *  SDEL 

S(2  = 

S(2  +  SG ( 2  j  *  SDEL 

p 

S(3)  = 

S ( 3)  +  SG(3)*SDEL 

C 

c 

STORE  THE 

SIGMA  MATRIX  FOR  1 

SK 

=  S(l) 

SDK 

=  S(2) 

SDDK 

=  S  (  3 ) 

c 

FIND  THE 

NEXT  VALUES  OF  THE 

c 

S (K+l) 

=  SPHI  *  S ( K) 

C 

c 

c 


220 


221 

222 

C 

C 


223 

C 

C 


c 

c 

c 


ZERO  A  TEMPORARY  MATRIX 

DO  220  L=1 , 3 

TEMPI (L)  =  0 
CONTINUE 
DO  222  L«l,3 

DO  221  M=1 , 3 

TEMPI (L)=TEMP1(L) 
CONTINUE 
CONTINUE 


+  SPHI ( L , M)  *  S (M) 


INPUT  BACK  INTO  SIGMA  MATRIX 
DO  223  N  =  1,3 
S  (N)  =  TEMPI (N ) 

CONTINUE 

STORE  THE  VALUE  OF  THE  S  MATRIX 

skpi  =  sm 

SDKPl  =  S ( 2 ) 

SDDKP1  =  S ( 3 ) 


NOW 


FIND 
PS  (3 
PS  (3 
PS  (3 
PS  (2 
PS  (2 
PS  l  2 
PS 
PS 


\  "  I  ' 


IF 


THE  F 

»  MATRIX 

)=PS  1 

;  3 , 3 ; 

\  - 

)  =PS  1 

;  3 , 2  ) 

)  ~ 

=  PSI 

;  3 .  i  i 

- 

)=PS  1 

!  2 , 3 ; 

— 

)=PS  1 

[2,2] 

— 

)  =PS  1 

1 2 , 1  ] 

— 

=  PS  I 

[1,3 

— 

)=PS  1 

1,2 

— 

)-ps(i,i: 

AT  STEP  K 


(9,  M 


'TIME  IS  ’ 
•A.B.C.D, ' 

A ,  B 
C,D 

RANGE  PHI 


450 


451 

C 


(HH  .LE.  0.0 \  THEN 
HH=50 
WRITE 
WRITE  ( 

WRITE  (9 
WRITE  (9 
WRITE  (9 
DO  450  L=1 
WRITE  (9 
CONTINUE 
WRITE  (9,*) 

DO  451  L=1 
WRITE  (9 
CONTINUE 
OUTPUT  THE  COVARIANCE 


PS(1, 
PS(1, 
PS(1, 
PS(1, 
PS  " 
PS 
PS 
PS 
PS  (1 


(1,3 

(1,2 


)  *  SG 

*  SG 

3) 

•SG 

3 

• SG 

2) 

)  *  SG 

2 

*  SG 

2) 

*SG 

1) 

*SG 

1) 

)  *  SG 

1) 

3 

*) 


,  TIME 


MATRIX' 


3 

M 


(RPHI (L,M) ,  M=1 ,3) 
SIGMA  PHI  MATRIX’ 
(SPHI (L , M) ,  M=1 , 3 ) 


MATRICES  AT  STEP  K 
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on 


401 

402 

300 

301 


WRITE  (9,*)  'THE  RANGE  COVARIANCE  MATRIX  IS:' 
DO  401  M=1 , 3 

WRITE  ( 9  ,  *  >  ( PR (M ,  N )  ,  N=l,3> 

CONTINUE 

WRITE  (9,*)  'THE  BEARING  COVARIANCE  MATRIX  IS 
DO  402  M=1 , 3 

WRITE  (9,*)  (PS (M ,  N )  ,  N-1,3) 

CONTINUE 
END  IF 

HH=HH-1 

RETURN 

END 

SUBROUTINE  2ERO(A,N) 

CLEAR  A  TEMPORARY  MATRIX 

REAL* 8  A (3 ,3) 

DO  301  L=1 , N 

DO  300  M-l.N 
A (L , M) =0 
CONTINUE 
CONTINUE 
RETURN 
END 
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C  THIRD  ORDER  MISSILE  SIMULATION  USING  CONSTANT  GAINS 

FOR 

C  THE  KALMAN  FILTER.  A  SECOND  ORDER  PROP  NAV  REFERENCE 

C  MODEL  IS  SIMULATED  FOR  PLOTTING  AND  PARAMETER 

COMPARISONS . 

INITIAL 

D  DIMENSION 

RNG (3) , S ( 3 ) , DELR ( 2 ) , GS ( 3 ) , GR ( 3 , 3 ) .RPHI (3. 3) ,SPHI(3,3) 

K=0 

NN=0 

C 

C 

c 

METHOD  RKSFX 

CONST  G=32 . 2 , D2R= . 0175, K2F=1. 66667 
TK  =  0.01 
MISSX0=0 . 0 
MISSYO=0 .0 
VM  =  2500.0 
AMO  =0.0 

READ  (2,10)  VT , AT , THDG , TGTX0 , TGTY0 
10  FORMAT  (F6.1,2X,F5.1,2X,F6.1,2 (2X.F10.2) ) 

TGTV=VT*K2F 
TGTA=~AT*G 
THDG=THDG*D2R 
TGTVXO=TGTV  *  S IN (THDG ) 

TGTVYO=TGTV*COS (THDG) 


C 

c 


INITIALIZE  THE 
RPHI (1,1) 


C 

C 


RANGE 
=  1 
TK 

.  5*TK*TK 
0 
1 

TK 
0 

......  0 

RPHI (3,3)  =  1 


PHI  MATRIX 


RPHI 

RPHI 

RPHI 

RPHI 

RPHI 

RPHI 

RPHI 


II  : 
i!  ‘ 


11: 
II. 
1!:?! 


INITIALIZE  THE 
SPHI (1 ,1) 
SPHI (1,  ‘ 
SPHI  1, 

SPHI (2, 

SPHI (2, 

SPHI (2, 

SPHI (3, 

SPHI (3, 

SPHI (3, 
CONSTANT 


GR (1.1) 
GR  < 1 , 2  ) 

811:11 

11:11 

3.1) 

3.2) 


2) 

3) 

1) 

2) 

3) 

1) 

2) 

3) 

STEADY 


BEARING  PHI  MATRIX 
1 

TK 

. 5*TK*TK 
0 
1 

TK 

0 

0 

1 

STATE  GAIN  VALUES  RANGE 


GR 
GR 
GR  ( 
GR  ( 


3  = 


C 

C 


C 

C 

ZERO 


GR (3,3)  =  0 

CONSTANT  STEADY 
GS  ( 1 )  =  1.5 
GS  2)  =  12.5 
GS  ( 3 )  =  1250, 


5 

0.0125 
0 
0 
1 
0 

0. 

24 


,0025 

.0 

1250 

.9991 


STATE  GAIN  VALUES  BEARING 


INITIALIZATION  OF  PROP  NAV  MISSILE  CONSTANT  VELOCITY, 
ACCEL 

LOS  =  AT AN 2 (TGTYO  -  MISSYO  ,  TGTXO  -  MISSXO) 

R= ( (TGTXO  -  MISSXO) *  *2  +  (TGTYO  -  MISSY0)**2  )**.5 
TTG0=  R/VM 
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PHDG* AT AN 2 (TGTYO+TGTVYO*TTGO-MISSYO , TGTXO+TGTVXO*TTGO-MISS- 
XO) 

VMXO  =  VM*COS (PHDG) 

VMYO  =  VM*SIN (PHDG) 

RKP1  —  R 

RDKP1  *  -VM*COS (PHDG-LOS)  +  TGTV*SIN (THDG) /COS (LOS ) 

RDDKP1  =  0 
SKP1  =  LOS 

SDKP1  =  (TGTVYO/COS (LOS)  -  VM* SIN ( PHDG-LOS )) /R 
SDDKP1  =  0 
BO  *  LOS 
BDO  =  0 
C 

RNG(l)  =  RKP1 
RNG ( 2 )  =  RDKP1 
RNG ( 3 )  *  RDDKP1 
S(l)  «  SKP1 
S  2)  =  SDKP1 

S  ( 3 )  =  SDDKP1 

C 

C  INITIALIZATION  OF  SECOND  ORDER  PROP  NAV  REFERENCE 

MODEL 

GO=LOS 
GDO  =  0 
C 

DERIVATIVE 

C 

C  TARGET  POSITION  UPDATING 

TGTHDG  =  ATAN2 (TVELX , TVELY ) 

TGTAX=TGTA *  COS ( TGTHDG ) 

TGT AY  = ( - 1 *  TGT A ) *SIN (TGTHDG) 

TVELX* INTGRL (TGTVXO , TGTAX) 

TVELY* INTGRL (TGTVYO , TGTAY ) 

XT=INTGRL (TGTXO , TVELX) 

YT=INTGRL (TGTYO , TVELY) 

C 

C 

C  THIRD  ORDER  PROP  NAV  MISSILE  POSITION  UPDATING 

BDOT  =  INTGRL ( EZ  C , BDDOT) 

B  =  INTGRL (BO, BDOT) 

AM  *  U 

MVELX  =  INTGRL (VMXO, -XM* SIN (PNHDG) ) 

MVELY  =  INTGRL (VMYO, AM* COS (PNHDG)  ) 

PNHDG  =  ATAN2 (MVELY , MVELX) 

PXM=INTGRL(MISSXO, MVELX) 

PYM=INTGRL (MISSYO , MVELY) 

C 

C  SECOND  ORDER  PROP  NAV  MISSILE 

GDDOT  *  -20*GDOT  +  ( SOLOS-GAMMA) *100 
GDOT  =  INTGRL (GDO, GDDOT) 

GAMMA  =  INTGRL (GO, GDOT) 

SOHDG  =  INTGRL (PHDG. 4* GDOT) 

SOXM  =  INTGRL (MISSXO , VM*COS (SOHDG) ) 

SOYM  =  INTGRL (MISSYO, VM*SIN (SOHDG) ) 

C 

DYNAMIC 

C 

C  THIRD  ORDER  PROP  NAV  MISSILE  GEOMETRY  UPDATE 

RM  =  { (XT-PXM ) *  *  2  +  (YT-PYM) **2) ** . 5 
LOS  =  ATAN2 (YT-PYM, XT-PXM) 

RDOTM  =  TGTV*SIN (TGTHDG) /COS (LOS) 

VM*COS ( PNHDG-LOS ) 

C  COMPUTE  THE  ERROR  TERMS 

DELR(l)  =  RM  -  RKP1 
DELR ( 2 )  =  RDOTM  -  RDKP1 
SDEL  =  LOS  -  SKP1 


CALL  KALMAN ( RNG , RM , RDOTM , RDDOTM , RK , RDK , RDDK , RKP 1 , 
+RDKP1 , RDDKP1 , DELR , S , LOS , LOSD , LOSDD , SK , SDK , 

+  SDDK , TIME , SKP1 , SDKP1 , SDDKP1 , SDEL , GS , GR , RPKI , SPHI ) 

C 

BDDOT=SDDKP1  +  10  * ( SDK-BDOT) +33 . 33333* (LOS-B) 

U  =  - 4  *  BDOT  *  RDOTM 
C 

C  SECOND  ORDER  PROP  NAV  MISSILE  GEOMETR.  UPDATE 

SOR  *=  (  (XT-SOXM)  **2  +  ( YT-SOYM)  *  *  2 )  *  *  .  5 
SOLOS  =  ATAN2 (YT-SOYM, XT-SOXM) 

SOU  =  — 4*GDOT*SORD 
C 

IF  ( PXM  .GT.  XT)  THEN 
CALL  ENDRUN 
END  IF 
C 

C  STATEMENTS  TO  SAVE  DATA  FOR  PLOTTING  WITH  DISSPLA 

IF  (K  .LE.  0)  THEN 

WRITE  (39,20)  PXM , PYM , SOXM , SOYM 
WRITE  (47,30)  TIME, LOS, U 
WRITE  (46,30)  TIME , SOLOS , SOU 
K=1G 

IF  (RM  .LT.  . 1 *R )  THEN 
K=3 
END  IF 

NN-NN+ 1 
ENDIF 
K=K-1 

20  FORMAT  (4 (2X.F10. 2)  ) 

30  FORMAT  (F5 . 2 . 2X.E11 . 3 , 2X, E14 . 6) 

SAMPLE 

C  STATEMENTS  TO  SAVE  DATA  FOR  GRAFAEL 

C  SAVE  (A)  0. I, XT, YT, PXM, PYM 

C  SAVE  (B)  0.1,  LOS , B , SK 

C  SAVE  (C)  0.1, LOSD, BDOT, SDK 

C  SAVE  (D)  0.1, SDEL 

C  SAVE  E)  0.1,RM,RK, RDDOTM, RDDK 
C  SAVE  ( F)  0.1,U,BDDOT 

C  SAVE  G)  0.1.GS1.GS2.GS3 

C  SAVE  (H)  0 . 1 , GR11 , GRl2 , GR21 , GR22 

C  PRINT 

0 . 1 , RM , RKP 1 , RK , RDOTM , RDKP1 , RDK , RDDOTM , RDDKP 1 , RDDK 
C  PRINT 

O.l, PNHDG , LOS , B , SK , LOSD , BDOT , SDK , LOSDD , BDDOT , SDDK 
CONTROL  FINTIM=1 0 . 0 , DELT= . 01 
C 

TERMINAL 

READ  (1,40)  MM 
WRITE  (1,40)  NK 
40  FORMAT  (F6.1) 

C  STATEMENTS  FOR  PLOTTING  USING  GRAFAEL 

C  GRAPH  (A/A,DE=TEK618)  XT 

( SC=1600 , LO=0 . 0 ) , YT ( SC=500 , P0=16000 ) 

C  GRAPH  (A/A,OV)  PXM  ( SC=1600 , AX=OMIT) , PYM ( SC=500 ) 

C  GRAPH  (B/B , DE=TEK618 )  TIME , LOS ( SC= . 025 ,LO=- . 1 > 

C  GRAPH  ( B/B , OV )  TIME (AX=OMIT) , B (PO=7 . 5, &C= . 025 , LO=- . 1 ) 

C  GRAPH  B/B,OV)  TIME ( AX=OMIT ) , SK ( AX=OMIT , SC= . 025 , LO=- . 1 ) 

C  GRAPH  (C/C, DE=TEK618>  TIME,l6sD 
C  GRAPH  C/C.OV)  TIME(AX=OMIT) ,BDOT(PO=7 .5) 

C  GRAPH  (C/C,OV)  TIME (AX=OMIT) , SDK (AX=OMIT) 

C  GRAPH  (D/D , DE=TEK618 )  TIME, SDEL 
C  GRAPH  (E/E , DE=TEK618 )  TIME , RM ( SC=2000 . 0 , LO=0 . 0 ) 

C  GRAPH  E/E , OV )  TIME  ( AX=OMIT) , RK ( SC=2000 . 0 . LO=0 . 0 ) 

C  GRAPH  ( F/E , DE=TEK61 8 )  TIME , RD60TM ( SC=500 . 0 ) 

C  GRAPH  ( F / E , OV )  TIME (AX=OMIT) , RDDK (AX=OMIT) 

C  GRAFH  (G/F , DE=TEK618 )  TIME,U 

C  GRAPH  H/G , DE=TEK618  TIME,GS1 

C  GRAFH  ( I/G , DE=TEK618 )  TIME , GS2 

C  GRAPH  J/G,DE=TEK618)  TIME,GS3 

C  GRAPH  ( K/H , DE=TEK61 8 )  TIME , GR11 
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TIME , GR12 
TIME , GR21 
TIME , GR22 
TIME , BDDOT 


C  GRAPH  (L/H, DE=TEK618 ) 

C  GRAPH  M/H,DE=TEK618 

C  GRAPH  (N/H , DE=TEK618 ) 

C  GRAPH  (O/F , DE=TEK618 ) 

END 
STOP 
FORTRAN 

SUBROUTINE 

KALMAN ( RNG , RM , RDOTM , RDDOTM , RK , RDK , RDDK , RKP 1 , RDKP1 , 

• 

RDDKP1 , DELR , S , LOS , LOSD , LOSDD , SK , SDK , SDDK , TIME , SKP1 , SDKP1 , 

*  SDDKP1 , SDEL , GS , GR , RPHI , SPHI ) 

C  SUBROUTINE  TO  ITERATE  A  KALMAN  FILTER  FOR  RANGE 

VARIABLES 

C  GIVEN  THE  COVARIANCE  MATRIX,  OBSERVATIONS 

REAL *8 

RNG ( 3 ) , RM , RDOTM , RDDOTM , RK , RDK , RDDK , RKP1 , RDKP1 , RDDKP1 , 

* 

DELR (2) , S ( 3 } , LOS, LOSD, RPHI (3,3) , A , B , C , D , TEMPI ( 3 ) ,SPHI(3,3) , 

*  GR (3,3) , GS ( 3 ) 

C 

C  FIND  THE  NEW  VALUES  OF  RPHI  FROM  THE  PREVIOUS  VALUES 

OF  SIGMA 


C 

c 


lie 

c 


120 

12- 


12: 

C 


C 

c 

c 


13  0 
131 

C 

c 


132 

C 


MATRIX 

DO  110  N  =  1.3 
TEMPI (N )  =  6.0 
CONTINUE 

CONSTANT  GAIN  INPUTS 
DO  121  L=1 , 3 
DO  120  M=1 , 2 

TEMPI (L)  =  TEMPI (L)  +  GR (L , M)  *  DELR (M ) 

CONTINUE 
CONTINUE 
DO  125  N=1 , 3 

RNG (N)  =  TEMPI (N)  +  RNG (N) 

CONTINUE 

SAVE  THE  VALUES  OF  RANGE  MATRIX  AT  STEP  K 
R0  =  RNG ( 1 ) 

RDK  =  RNG (2) 

F1DK  =  RNG (3) 

ZERO  THE  OLD  RANGE  TEMPORARY  MATRIX 

FIND  THE  ESTIMATE  OF  THE  STEP  K+l  FOR  THE  RANGE  MATRIX 


DO  131  L  =  1,3 
D-  130  M  =  1,3 

TEMPI (L)  =  TEMPI (L) 
CONTINUE 
CONTINUE 


+  RPHI (L,M)  *  RNG (M) 


SAVE  THE  VALUES  OF  RNG (K+l/K) 
DO  132  N=1 , 3 


RNG (N) 
CONTINUE 

RKP1 


TEMPI (N) 


RNG ( 1 ) 

R0KP1  =  RNG  2 
RLDKP1  =  RNG (3) 

TO  ITERATE 


A  KALMAN  FILTER  FOR  SIGMA 


C 

C  SUBROUTINE 

VARIABLES 

C  GIVE:;  THE  COVARIANCE  MATRIX,  OBSERVATION. 

C  WHERE  H=  (1  0  0) 

C 

C  NOW  FIND  THE  CURRENT  VALUES  OF  THE  SIGMA  MATRIX 

S  a)  =  S  (1)  +  GS  (1)  *  SDEL 
S{2)  =  S  2)  +  GS ( 2 ) *  SDEL 
St3)  «  S(3)  +  ' 

C 


GS (3) *  SDEL 
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C  STORE  THE  SIGMA  MATRIX  FOR  USE  IN  THE  PROGRAM 

SK  =  S ( 1 ) 

SDK  =  S  2 
SDDK  =  S ( 3 } 

C  FIND  THE  NEXT  VALUES  OF  THE  SIGMA  MATRIX 

C 

C  S ( K+l )  *  SPHI  *  S (K) 

C 

C  ZERO  A  TEMPORARY  MATRIX 

C 

DO  140  L=1 , 3 

TEMPI (L)  =  0 

140  CONTINUE 

DO  142  L«l,3 

DO  141  M-1,3 

TEMPI (L)=TEMP1 (L)  +  SPHI(L,M)  *  S(M) 

141  CONTINUE 

142  CONTINUE 
C 

C  INPUT  BACK  INTO  SIGMA  MATRIX 

DO  143  N  =  1.3 
S (N)  =  TEMPI (N) 

143  CONTINUE 
C 

C  STORE  THE  VALUE  OF  THE  S  MATRIX 

SKP1  =  S (1 

2 
3 
C 

RETURN 
END 

SUBROUTINE  ZERO (A, N) 

C  CLEAR  A  TEMPORARY  MATRIX 

C 

REAL*  8  A (3, 3) 

DO  201  L=1 , N 

DO  200  M«1,N 
A (L,M) =0 

200  CONTINUE 

201  CONTINUE 
RETURN 
END 


SDK? 1  =  St 

SDDKP1  *  S( 


) 
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